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MODELS FOR FORCE MOTION AND 
INTERACTION THAT REPRESENT 
UNCERTAIN PERCEPTION 


D. P. Gaver 
P. A. Jacobs 
M. A. Youngren 
S. H. Parry 


SUMMARY 

This report records work carried out on a theater-level simulation modeling effort 
entitled J-STOCHWARS (JS). It is a connected set of working papers: minimal attempt 
has been made to polish the prose, and redundancies occur. There remain needed 
additions and modifications. These will be made as the authors adapt the ideas for use in 
planning models that are under current, and future military development. Work on the 
simulation model has been supplanted by an analytical scoping model, Battlespace 
Information War (BAT/IW). 

The emphasis of the report is upon recognition, characterization and exploitation and 
control of aspects of uncertainty in military operations. It emphasizes simple low- 
resolution representations of the products of generic imaging sensors. Later work will 
deal with other input, such as ELINT and COMINT. We provide ideas as to how such 
information products can be combined or fused, basically using a Bayesian and 
Gaussian/Normal-distribution theory format for computational and conceptual 
convenience. We investigate processes for probabilistically updating opponent course of 
action (COA) estimates. 

Many of the computational ideas and algorithms have been incorporated in a 
computer program termed Joint Warfare Analysis Experimental Prototype (JWAEP). The 
latter runs on a SUN computer. 
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1. Background 

It is a truism to state that when several military forces jointly interact within a theater 
to oppose the actions of another they do so under conditions of substantial uncertainty as 
to the outcome. In fact, Carl von Clausewitz has said “War is the realm of uncertainty: 
three quarters of the factors on which action in war is based are wrapped in a fog of 
greater or lesser uncertainty. A sensitivity and discriminating judgment is called for; a 
skilled intelligence to scent out the truth.” The purpose of this report is to list many of the 
sources of that uncertainty and to explore ways of describing their possible impact on a 
theater-level conflict, as the latter may develop over time. The style used is that of 
quantitative mathematical modeling, but many of the numbers {parameters ) are at best 
tentative, as are the ways of representing their interaction (the models). Mathematical 
modeling requires specificity to obtain actual numerical “results”, but in many situations 
these results must be treated as suggestive and tentative, and neither accepted nor 
discounted too hastily, without assessing both vices and virtues of the modeling option 
currently in use. Our approach will be to propose and give arguments for, and against, a 
variety of such options for many aspects of theater-level conflict problems without 
exclusively endorsing any one. 

Our purpose is to develop models that have the “feel” of aspects of real campaigns, 
but that omit many details. They are “low resolution” by design. In particular, the optional 
approaches should allow the user-analyst to efficiently experience alternative futures that 
might plausibly follow from adopting certain doctrinal options. 

This report is a record of work done on various aspects of theater level combat 
modeling; the emphasis is on simple perception modeling and not on attrition as this 
might be linked to perception. Work on the latter is reported elsewhere. 
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1.1 Sources of Uncertainty 

Sources of uncertainty in military operations, and options and variations in the 
modeling thereof, are many. For a recent general review see Hughes (1994). Prominent 
among them is the actual, operationally relevant, status of various physical aspects of the 
theater environment, e.g. of the geographical terrain and weather, as influenced by season 
and time-of-day; others can be identified. Abbreviate designation of this class by PEUF 
for physical environment uncertainty factors. PEUF influence speed and predictability of 
force movement, visibility and hence detection and classification of force elements, 
fortification potential or the defensibility of specific locations, communications, 
equipment performance and reliability, ultimately even morale, all to varying degrees, 
many of which may be anticipated by an opponent, but with inevitable error. 

Another prominent and related uncertainty source has to do with the trustworthiness 
of intelligence, this term generically referring minimally to knowledge of the present 
location, activity and current motion of both opponents’ forces by their respective 
opponents. This knowledge, with its inevitable uncertainty, depends on the types, manner 
of utilization, and product of information-gathering sensory assets. These range from 
human scouts and informants (HUMINT) to electronic intelligence (ELINT) acquisition 
assets such as radars or sonars, either stationary or carried on various mobile platforms 
(manned or unmanned reconnaissance, aircraft or ground vehicles, satellites, etc.) that 
actively search certain aspects of terrain for opponent asset occupancy, but with 
unavoidable omissions and inaccuracies. 

Additional sources of information, also susceptible to errors and omissions of 
acquisition and interpretation, including active opponent deception, are a variety of assets 
such as satellites, manned reconnaissance platforms, UAVs, and TUGVs. All of these 
ultimately supply images that must be interpreted. We group these in the category 
abbreviated as MINT. In addition, there is the information to be gained from sensing the 
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communication emissions of an opponent force and interpreting their meaning. This area 
is called COMINT, and is useful but also susceptible to errors of misinterpretation and 
delay. The uncertainty associated with (only!) the physical performance of information¬ 
gathering assets, e.g. its probability of detection of a specified target asset such as a tank 
(or tank company) or armored personnel carrier, its probability of correct classification 
given detection, its accuracy of location of the sensed asset (conditional on the sensor 
system operator’s skill and the surrounding physical and operational conditions) can be 
abbreviated as PIAU, or physical intelligence asset uncertainty factors. 

When actual combat occurs further sources of uncertainty must be confronted, namely 
the physical capabilities of the weapons brought to bear. Even if intended targets are 
well-located and weaponry is extremely accurate (“smart”) and well-matched to target 
vulnerability there are ample opportunities for errors and equipment (and operator) 
failures that can lead to partially- successful mission accomplishment, and hence 
requirement for follow-up action. We lump such physical weapon asset uncertainty into a 
category abbreviated by PWAU. It is apparent that even smart weapons can be decoyed or 
otherwise duped, which adds a further element of uncertainty, and necessitates thoughtful 
battle damage assessment (BDA), in turn influenced by PIAU as above. 

The above rough categorization of the major physical components of uncertainty that 
must be considered by both opponents and allies in theater level warfare is not 
exhaustive. There are other far more tenuous components that are, however, significantly 
affected by the stark physical influences mentioned above. A major component of the 
overall uncertainty must be that concerning the territorial objectives of, say, an aggressor 
into a neighboring area, and the willingness or resolve to attain those objectives at the 
expense of losses of assets and his/her own territory. It is far more difficult to credibly 
quantify the uncertain tenacity of an enemy force in the face of actual or prospective 
attrition than it is to model physically-based uncertainties associated with unopposed 
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enemy advance rate through somewhat unknown, and unexpectedly variable, terrain 
while accounting somewhat approximately for that force’s size and weapons and 
intelligence capabilities. Opponent resolve and tenacity have often been modeled by 
thresholds that limit attrition in case combat occurs, perhaps by withdrawing, but also by 
calling for reinforcements, which are then unavailable to another area or conflict in a 
different sub-theater. Note that the thresholds can be of various forms: they can reflect 
own force losses (or prospective losses), or, additionally, estimates of opponent force 
losses; the latter may well be “known” only roughly. Thresholds may also be expressed as 
current or projected rates or estimated rates (derivatives) of loss, cf. Helmbold (1971). 

A generic scenario of Red and Blue forces moving towards a common objective in an 
invasion, occupation, protection scenario is described in the next section. A brief 
description of the physical theater in terms of physical nodes and arcs is given. 
Increasingly detailed models of force movement along arcs are presented in Sections 2 
and 3. Models for detection and tracking of units are also introduced. Procedures for 
estimating each side’s units’ velocities, arrival time at the objective and unit size are 
proposed and studied. In Sections 4-6, procedures for estimating the number of assets 
on a node or arc based on binomial, beta-binomial and multinomial models for sensor 
observations are proposed and studied. In Section 7, procedures for estimating the 
number and types of military units on a node/arc based on estimates of the number of 
assets at that node and neighboring nodes are studied. Section 8 describes a formulation 
of each side’s course of action (COA) in a conflict. Procedures for one side to estimate, 
using sensor observations, the probability the other side is following each COA are 
described. 

The present report is a progress statement that documents ideas and concepts that may 
be useful in other modeling projects. Some, but not yet all, of these ideas have been 
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realized in JWAEP (Joint Warfare Analysis Experimental Prototype), a demonstration 
software package described in the JWAEP User’s Manual; cf. Youngren (1996). 

2. Invasion-Occupation-Protection Dynamics: A Simple Prototypical 
Example 

We propose for initial discussion a simplified prototype for many real situations: 
Iraq’s invasion of, and expulsion from, Kuwait, or its possible invasion of other 
neighboring countries; North Korea’s possible invasion of South Korea, and perhaps 
others on a smaller scale. It is stipulated to possess these elements: 

(a) A region , abbreviated % that is in dispute, possibly because of its desirable 
resources (oil, resort areas and gambling casinos, bountiful production of, say, hog 
jowls or titanium). 

(b) One (or more) unfriendly invading forces, termed Red, that intend to capture 
territory and exploit resources of % They must enter the region by crossing its 

boundaries. Some entry points have advantages, i.e. are less well-protected or 
defensible than others. Let R(t) be the size of Red force(s) at t. R(t) is actually a 
vector whose components are numbers of force types and their locations; R(t) is Red 
ground truth at time t. 

(c) A protective or defensive force, or forces, that also originates outside of W; call it 
Blue ( B ). Its force size is B(t). If an offensive move by R occurs, B responds in one 
of several ways. B(t) is also a vector as is R(t); it is Blue ground truth at t. There is 
dynamic interaction between R{t) and B(t). 

(d) A civilian or non-combatant indigenous population', it may well be that the 
population must be split into segments sympathetic to R, C/j(t) in number, and those 
sympathetic to B, Cs(t), with total C(t) = C/j(t) + Cg(r). Finer distinctions are 
always possible, in which case, the above state variable may be a vector. In many 
situations affiliations would change as combat progresses and population elements 
alter their perceptions as to whom they might wish to ally themselves with. 
Considerations for such a population can well be a major concern in operations- 
other-than-war (OOTW) or peacekeeping scenarios. 
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(e) An indigenous military force able to delay, but not repel, an unfriendly invading 
force, or forces. Let I(t) generally denote the size of that force at time t; realistically 
it will be vector-valued so as to represent force components: infantry, armor, 
artillery, air, etc. There may also be indigenous forces that are allies of the invaders; 
let the number of these be Ig(t). Those allied with B are I sit) in number. 

Of course the implications of the above terminology need not hold for all situations, 
i.e. that the region, % is in friendly (to B ) hands initially and is threatened by “bad guys” 
and protected by “good guys”. In fact, W may be initially in the hands of R, and its 
domination could come under attack by B. And the geography must be specified in 
enough detail to portray a region and the relevant surrounding territory. If desired, regions 
can be represented as a union of non-overlapping subregions. 

Geography: Arcs and Nodes 

A theater, that is the geographical entities relevant to the elements (a) - (e) above, is 
conveniently represented by a system of nodes : physical locations, fixed in space, and of 
importance either because they have symbolic value, e.g. are centers of governmental 
authority, or are suitable for defense by fortification or asset concealment because of 
topography, or are significant locations for military assets such as airfields, surface-to-air 
missile sites, (relatively) immobile command centers or possibly C3/I facilities. Nodes are 
viewed as interconnected by arcs : road systems or other routes for the transit of ground 
combatants. For some purposes it is useful to represent arcs as a sequence of discrete 
locations (Subarcs, or compartments) between which military units transition. 

In a naval context it may often be convenient to specify geographical locations as 
being within particular broadened arcs or spatial compartments, shaped as squares or 
hexagons. 

A primitive arc-node setup appears below: 
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Arc 1 (Red’s path to 5V) 


Arc 2 (Blue’s path to 5V) 



Node 1 (Red Force; Node 2 Node 3 (Blue Force; 

Home Location) ^V; (Disputed Region; Home Location) 

Civilian & Indigenous 
Force Location) 

Figure 2.1 

The above simplistic setup can be elaborated greatly: realistically there may be many 
routes from Node 1 to Node 2, and Node 3 to Node 2, possibly by way of intermediate 
nodes. 

The arc-node representation is perhaps more appealing when alternative routes can be 
cleanly and crisply distinguished because of topography (a substantial ground force does 
not easily traverse mountain ranges) than when the topography is bland and neutral, as in 
a desert, or in the open ocean. From the point of view of air warfare, wherein aircraft 
from a source node may transit to a target area to conduct reconnaissance or an offensive 
strike, arcs are a less attractive representation than is a broad “tiling” of the region, 
possibly by square or hexagonal subregions or compartments. Sorties of aircraft simply 
move, in appropriate time, from tile or compartment (center) to a contiguous tile or 
compartment (center) as they transit a region. Alternatively, straight line flight paths 
could be scripted through the region. 

Dynamics 

Military operations by, say. Red, and in response by Blue, can be specified in terms of 
arcs and nodes. For example, in Figure 2.1 a Red force moves on Arc 1 from N1 to N2. 
Some options for the dynamics are listed below. 
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Option 1: Simple Random Time-On-Arc (Time-On-Node) 

It is a reality that even an unopposed ground force unit will move across a given 
terrain route in a somewhat unpredictable, variable, or uncertain time. This is one 
consequence of PEUF (Section 1). The simplest way of representing one such transit time 
is to let it be a single random draw from an appropriate population, or realization, T, of a 
random variable, T. Then if D is the distance, say from Node 1 to Node 2 in Figure 2.1, 
the velocity of transit is V- D/T, a realization of a random variable. Having obtained Tit 
is then easy to locate the unit’s location (center) at time 0 <t<T after it leaves the origin 
node if we assume a constant velocity: its distance from the starting node in the direction 
of motion is D(t) = Dt!T\ this number locates the unit on the arc and is referred to as 
ground truth. Note that actual ground truth will vary across realizations of the model 
because of the random variability of T. Of course the distributional properties of T will 
potentially depend on terrain type, even direction of motion, (e.g. up or down hill), unit 
type, time of day, season of year, and perhaps other factors as well. Convenient candidate 
distributional forms are the Gaussian or normal (adjusted to truncate away possible 
negative values), the inverse Gaussian (which has naturally positive support), the 
lognormal, log-logistic, gamma, Weibull, or other distributions with positive support. A 
practical issue is that of actually specifying numerical parameters to specify the 
distributional form selected. Estimates of the mean and variance should be available from 
field data, but analyst or other expert judgment may well be required, at least initially. 
The effect of concomitant explanatory variables can be incorporated by regression 
modeling; see McCullagh and Nelder (1993). It will thus be possible to represent weather 
effects by regression. 

The advantage of this approach is its simplicity. Disadvantages or unrealities include 
the fact that a unit’s velocity across an arc is represented as constant once the time T is 
selected; of course this will change when another plausible history is selected, i.e. on 
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another model evolution. If an event occurs while the unit is in transit, e.g. if it encounters 
an obstacle such as a mine field and waits to sweep it, or alternatively bypasses it at the 
cost of additional delay, a new random number must be drawn. Alternatively, the mean 
and variance of the original time, T, may be adjusted to reflect such extra delays. Also, if 
the unit in transit encounters an opponent and pauses for combat, another random draw 
will be required to represent the future motion of the unit. Note that for the purpose of 
reducing computing demand all elements or Subunits (e.g. companies) of a unit (e.g. 
brigade) are represented as moving together. They may be conveniently viewed as located 
in an assigned formation (template) within a moving rectangular matrix. 

Option 2: Markovian Movement Through Subarcs 

To better represent potential local in-transit variability of motion, suppose that an arc 
between two consecutive nodes is divided into several consecutive Subarcs or 
compartments; see below, where there are exactly k Subarcs. 



Figure 2.2 

Imagine that a Unit (e.g. brigade) or collection of Subunits (e.g. companies composing a 
brigade) leave Ni and enter SAi at t = 0; they remain in SAj for a random residence time, 
T\, and thence move to SA 2 , from which they transit to SA 3 after random residence time 
T 2 , and so on through SA* and to N 2 . If the above process is defined in continuous time 
(t e R + ) and the transit/ residence times are independently and exponentially distributed 
with parameters Xj then the location of a unit, i.e. the Subarc occupied, at time t is a 
Markov process in continuous time. If the process is defined in discrete time, and the 
transit/residence times are independently and geometrically distributed with probability of 
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one-stage advance being pj (and staying in place on SAj with probability 1 - pj) then the 
location of the unit at time (t e (0,1,2,...) is Markov in discrete time. 

Example: Suppose the motion process is in continuous time and all Subarc times, T\, 
Tj, ... are independently and identically exponentially distributed, with parameter X. 
Then the unit moving is located at/within Subarc j (j < k) at time t with the Poisson 
probability 


.P{Unit in subarc j at t] = e h 



( 2 . 1 ) 


the most likely Subarc to be in at t is j = [Xt\\ if [Xt\ >k + 1 the Unit has arrived at N 2 . 
Now to calibrate parameters X and A: to a specified mean unit velocity V at which the 
distance D is covered we require 


k_ D 
X~ V 


T. 


( 2 . 2 ) 


Suppose in addition we want to specify that the standard deviation of T is a fraction,/ of 
T . This means that 


= (2.3) 

which implies that k= l// 2 . Thus if /=0.10 the number of Subarcs k= 100, while if 
/= 0.25, k= 16: the larger the variability the smaller the number of Subarcs required. 
With a moderate number of Subarcs (even 16) it will usually be adequate to use the 
normal approximation to the Poisson. The parameter X, or rate at which a Subunit moves 
off a Subarc is, from (2.2), obtained as X = k/ T , or, more intuitively, MX = T Ik. 

An extension of the previous simple model is to allow backward as well as forward 
motion, and also pauses on an arc of more than one basic exponential duration. 


11 



Numerical Illustration 2.1. Independent Simultaneous Approach; Race Towards a 
Common Objective. 

Figures 1-4 show the results of a simulation of the simultaneous approaches of Blue 
and Red forces towards the common objective J\ r = N 2 . For this illustration the two forces 
move independently of each other. The Blue force starts from a distance of Dg = 200 
miles and moves at mean velocity Vg = 15 mph; Red from Dr = 100 miles, mean velocity 
10 mph. Each distance is divided into 20 Subarcs. There were 1000 replications. A force 
is said to be late if it is not the first to arrive at % The resulting lateness time is the length 
of time after the first arrival until arrival of the late force. Of special tactical/operational 
interest is the lateness of either Blue or Red: for the present situation Blue’s lateness 
probability is estimated to be about 0.80; if late, Blue’s lateness time mean is 4.4 hours, 
but the distribution of lateness time, given Blue arrives last, is substantially spread out, 
providing Red considerable time to prepare a position on N 2 . On the other hand, Red is 
late with probability about 0.20; if late, Red’s lateness time mean is only about 2 hours 
and seldom more than 5 hours. Note that a deterministic assessment of arrival order 
would be based on comparison of Db/Vb = T b = 200/15= 13.3 and Dr/V r = T r = 

100/10= 10, so on that basis Blue would always arrive 3.3 hours after Red. A more 
refined and easy calculation is based on a normal approximation 




Discussion. The above illustration is quite rudimentary and unrealistic in that no sensor 
activities and resulting information, and subsequent action changes, are modeled. It does, 
however, illustrate the effect of variability in transit. 

Numerical Illustration 2.2. Motion with Detection and Estimation of Arrival Time. 

Consider this movement scenario, similar to the above, but modeling detection by 
Blue. Once again, Red is initially located a distance Dr from the objective % Blue is 
initially located a distance Dr from the objective % the mean velocity for Red 
(respectively Blue) is Vr (respectively V B ). 

Assume the distance Dr (respectively Dr) is divided into Nr (respectively Nr) equal 
Subarcs. Assume the time to transverse Subarc i for Red (respectively Blue) Tr(i) 
(respectively Tr(i)) has a distribution with mean (Dr/Vr )/Nr (respectively ( Dr/Vr )/Nr). 
Example. Dr = 100 miles Dr = 200 miles 

V R =10 miles/hr. Vr =15 mph 

Both Blue and Red traverse 30 Subarcs. Suppose now that Blue can detect Red. The 
probability that he detects Red on a Subarc, given her presence there, is 8, temporarily 
assumed constant. Thus it may be representative of wide-area theater-level, but not unit- 
level, overhead assets, implicitly cueing a “soda straw” sensor such as UAV. The number 
of the Subarc on which first detection is made has a geometric distribution (if there were 
an unlimited number of subarcs) 


P{M = k} = S(l-5) k ~ l k = 1,2,.... 

Assume that Blue tracks Red for d t Subarcs after detection at M. Blue’s estimate of the 
velocity of Red is 


Vr = 


M+d, 

2X0 

1= M±\ 

d,D R /N R 


-l-i 


( 2 . 6 ) 
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where Tr(i) is a time to transit Subarc i, a realization of Tr(i). Blue’s estimate of Red’s 
time of arrival at 5Vis 


Ar = Dr 


1 - 


M + d t 
Nr ; 


a ‘ V R ,=M+1 


(2.7) 


M+d, 


The elapsed time since the start of the replication is T R (i ). The number of the 

/=i 


Subarc that Blue is on at this time is 


( M+d, k 1 

k: ( 2 . 8 ) 

where Tg(i) is Blue’s transit time on his Subarc i. Given Blue’s estimate of his own 
velocity, Vg, Blue can estimate an approximate probability that Red will arrive at (N x 
time units before Blue will. 


1~Pl{x) 


P\ 


N B 

^Tb(})> Ar + 


[i=N E (B )+1 



(2.9) 


If the Tr( i) are exponentially distributed then the above calculation involves the survivor 
function of a gamma random variable with shape parameter Ng - Ne{B). 

Blue can potentially use the above results to calculate the chance of successfully 
carrying out certain actions to slow or damage the Red force, e.g. before that force 
reaches 5V; Blue may speed up, send strike or assign indirect fire assets to attrit and delay 

Red, or use assets such as FASCAM (Family of Scattered Mines; an air or artillery- 
delivered minefield) to create obstacles in Red’s path. 

Numerical Illustration 2.3. Example. 

Figures 5-10 display results for a simulation with 1000 replications in which the 
time to transit a Subarc is normally distributed. The number of Subarcs is 20 for both Red 
and Blue. The mean time for Red to transit a Subarc is (Dr/ Vr )/20 with standard 
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deviation (0.3 ){Dr/V r )/20. The mean time for Blue to transit a Subarc is (Dg/P# )/20 with 
standard deviation (0.3)(Z)g/Fj})/20. Figures 5 and (respectively 6) show histograms of 
the amount of time to transit the arc for Red (respectively Blue). 

The probability of Blue detecting Red on a Subarc is 5= 0.3. Figure 7 displays the 
number of the Subarc on which Red is detected. 

Blue then uses two additional Subarcs to estimate the velocity of Red. There are 3 
replications in which Red is first detected on Subarcs 18-20. For these replications Blue 
does not estimate the velocity of Red. Figure 8 displays a histogram of Blue’s estimates 
of the velocity of Red. Figure 9 displays a histogram of Blue’s estimate of the mean time 
remaining until Red achieves his objective. Figure 10 displays a histogram of Blue’s 
estimate of the probability he achieves the objective first. 

Numerical Illustration 2.4. Velocity Estimation by Exponential Smoothing. 

Figures 11-14 show the results of a simulation in which there are 1000 replications. 
In each replication Red is Dr - 100 mi. away from his objective and is traveling with a 
mean speed V R = 10 mph. The arc along which Red is traveling is evenly divided into 20 

Subarcs; thus each Subarc has a distance Dr(S) = 100/20 = 5. The time Red takes to 
travel a Subarc is normal with mean D R /{Vr x 20) = 0.5 hours and standard deviation 

0.2 x 0.5 hours; the normal random numbers are left censored at 0.05. The probability 
Blue detects Red on a Subarc is 0.3. Once Blue detects Red he is able to track Red. 

One procedure for Blue to use to estimate Red’s velocity once Red is detected is 
exponential smoothing. Let V(i) be the estimate of Red’s velocity after Red passes 

through Subarc i after Blue detects Red. Blue’s estimate of Red’s velocity after Subarc 
i +1 is 


V{i +1) = (l - a)V(i ) + a 


M[+i) 

Ms) 


15 



where 7/j(* + 1) is Red’s time to traverse Subarc i + 1. In the simulations a = 1/3. The 
estimated velocity can be used to estimate Red’s time of arrival at the objective. 

Figure 11 displays the histograms of the estimate of Red’s velocity. One histogram 
displays the estimate obtained by dividing the distance of two Subarcs by the time Red 
takes to traverse the two Subarcs after the Subarc on which he is detected. The other 
histogram displays the histogram of the exponential smooth estimate using those Subarcs 
after detection that were traversed before time 6; the time of 6 was chosen arbitrarily as a 
perception update time. Those replications for which there are not two Subarcs after 
detection available for estimation before time 6 are not used. As expected, the 
exponential smoothed estimate has less variability than the two Subarc estimates. 

Figure 12 displays histograms of the error of the estimate of Red’s time of arrival at 
his objective. The estimate using the exponential smoothed estimate of velocity is less 
variable than that for the velocity estimate using two Subarcs; this is to be expected. 

Figures 13-14 display results using the exponential smooth estimator of velocity at 
times 3 and 6. In each case a replication is deleted if it does not have two Subarcs to use 
before reporting. Figure 13 shows the histograms of the velocity estimates and Figure 14 
shows histograms of the errors of the estimates of Red’s time of arrival at the objective. 
The additional observation time to estimate the velocity of Red not surprisingly decreases 
the variability of the estimates. 

Example 2.5: Unit Markov Motion 

Suppose, for specificity, that a Red Unit of population size R (measured in number of 
Red Subunits, e.g. companies if R is one or more regiments or brigades) enters Arc 1 
between Ni and N 2 (Figure 2.2) via SAi at t = 0: Let i?i(0) = R and let Rj(t) = Rii the 
Unit is on Subarc j at time t; otherwise Rj(t) = 0. Represent the Unit’s motion as 
probabilistic as follows: 
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PjRed Unit moves from SA y - to SA y+ i in time step [t,t +dt)\x R (t) = xj 
= rj(t,dt;x) 

otherwise it remains at SAj for (t, t + dt) with probability 1 - rft, dt; x). The time 
increment dt need not be infinitesimal, e.g. it may be one hour. 

The (vector-valued) variable Xft) represents a general external 
influence, based on perception, upon the probability of advance by 
Red. There exists a corresponding external influence variable, Xs(t), 
for Blue. 

For instance Xgtf) may reflect the perception of Red that a Blue force of some 
estimated size is concentrated at a corresponding Subarc of Arc 2 (Figure 2.1) and that 
the corresponding probability of advance to disputed Node 2 has a specified estimated 
current value; this might be such that the Blue force will reach N 2 (-3V) within a specified 

(small) number of time steps with high probability. If this results in Blue reaching N 2 
sufficiently in advance of Red to fortify and hence obtain an advantage. Red must decide 
among operational options, such as (a) continue as planned; (b) speed up and reach N 2 
first (with high probability) and itself fortify N 2 ; (c) place obstacles in Blue’s path (e.g. 
emplace a minefield, or setup an ambush, at a Subarc between Blue’s current location and 
N 2 so as to delay Blue’s progress, simultaneously proceeding towards N 2 ; (d) other (many 
possibilities!). Blue has comparable operational options. 

Example 2.6: Suppose there are Cft) Subunits on/in Subarc j at time t; j <k in 
Figure 2.2. Then a number, Eft, t + dt), may be ordered to move (emigrate) from Subarc j 
to Subarc j + 1 in time dt (for the present dt is a fixed non-zero time interval analogous, 
but not necessarily equal, to the PRT (perception review time) A). Likewise, a number 
Ej.\(t, t + dt) on j- 1 may emigrate from Subarc j- 1 to Subarc j, leaving a net number on 
Subarc j equal to Cft + dt) = Cft) + Ej.\(t, t + dt)- Eft, t + dt). The above transfer of 
Subunits is only feasible if there is a positive number of Units or Subunits on Subarcs j- 1 
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and j to immigrate and emigrate, which means that both Ej.\ and Ej are functions of their 
parent Subarc population levels, Cj.\(t) and Cj(t), being zero if those population levels are 
zero. 

Example (Continued): Suppose that for theSubarc 


Ej(t,t + dt) = Cj(t)(l - e ~^ dt ); 1 < j < k. 


( 2 . 10 ) 


Then 


Cj(t +*) = Cj(t) + C ; -,(/)(l - A"")- Cj(t)[ 1 - e-W) 
= Cj(t)e-^ + C H (t){ 


( 2 . 11 ) 


which is feasible, giving a positive result. Furthermore it is readily advanced in time, e.g. 
starting from the initial condition Ci(0) = C, where C is the population size of a Unit 
about to move onto the arc at t = 0. 


C k ( t + 2 dt) = C k ( t + dt)e~^ dt + Ck-\ {t + dt){ 1 - e 4k - ld ‘ ) 

= [C* (*>-*** +C k .,{t){\-e-^ dt )Y^ kdt 

+[c k -,{t)e-* k -' dt + C k . 2 {t)( \-e~^ dt )](l-e-^ dl ) (2.12) 

= C k (t)e~ 24kdl + Ck-i (/)[(l - e^ dt )e~ 4kdt + (l - e~^ dt )e~ 4k - ]dt ] 

+C k - 2 (0(l - e~* k ~ 2d ‘ )(l ~ e ~^ k -' dt ). 

Example (Continued): If = £ i.e. if the transfer probability is the same for all Subarcs 
it can be seen that, putting ndt -1, 


CK0-C t (**) = tc y (0)f " ( 2 . 13 ) 

7=0 V* J J 


In words, if there are initially C/0) Subunits on/in Subarc j (0 <j< k) then any one of 
these must make exactly k-j transitions between intervening Subarcs in n time steps in 


18 


order to be on/in Subarc k>j at time ndt\ there are n ways of selecting those that 

“make it”. The formula (2.7) shows that if Co(0) = C Subunits start at time 0 then the 
number at exactly the Subarc at time ndt = tis 

= (2.14) 

The Subarc containing the maximum number of Subunits is approximately k = «(1 - 
e -< 5 *) when n is large; the number of Subunits on that particular Subarc approaches zero 
as n becomes large, at a rate proportional to l/Jn, and those nearby have about the same 

number: Subunits become quite dispersed. 

Comment. The above transition law, with £ or constant, could be appropriate for a 
disorganized force advance — or perhaps better, retreat — but is not truly representative 
of the motion of a structured force’s coordinated advance towards a territorial objective. 
However, the dispersive effect noted above can be remedied by introducing appropriate 
time dependence into the parameters 

3. A Red-Transit, Blue-Perception Vignette 

At time t = 0, a Red Unit (e.g. Army brigade), composed of C Subunits, moves off a 
physical node onto a transit node or arc. It spends a random transit time, T, on the arc, 
after which it reaches a subsequent physical node. Blue’s objective is to estimate current 
force size, C, on the arc from observations available from a sensor system. 

Opponent (Blue) Perception: Observational Modeling 

Suppose a Blue Overhead Sensor Suite is constantly observing the arc on which the 
Red Subunit is moving, and that the random time, D, to detect any one of the Red Unit’s 
Subunit components that leaves in a Lost or untracked state has exponential distribution 
with mean 1 IS (rate 8). Having detected such it tracks that Subunit for a random time L 
until loss of track occurs; L is assumed exponential with mean I/ 77 . Search resumes, and 
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the unit is detected again at a time that is an independent replica of D (presuming it is still 
on the arc); it is again susceptible to loss, and so on. The alternating sequences {D^} and 
{Lk} are assumed independent. Furthermore, individual Subunits are assumed to be 
detected independently. It is clear that the appearance on the arc of the C Subunits 
triggers C independent 2-state Markov chains. Note that a Subunit can be in a detection 
state at t = 0, as it begins transit, although the implication of the above discussion is that 
some Subunits begin their transits anonymously, i.e. in a Lost or undetected state. It 
follows that the number of units that start transit in a Lost state, C/,(0), and those that start 
in a Detected/Tracked state, Q)(0), are, at Blue’s (first) Perception Review Time (PRT) 
t = A, in either a Lost or Detected/Tracked state Cld( A), Cll( A), Cdd( A), Cdl( A), 
obeying the dynamics of a simple Markov chain with transition probabilities 

^{C iD (A)= v„(A)} =( Cl(0) ]p ID (A)''“ {4) (l-to(A)) C ‘ (0) "'“ (4) (3-D 

IM A )J 

/>{C It (A)= v u (A)} =f Cl<0) L(Ap (4) (l-MA)) Q<0, ‘ , ' 1 ‘ ,4) (3.2) 

v Vzx (a)> 

P{C DD (A) = Vdd(A)} = [ ^ L c (A)"“> (i >(l- p DD (t,)) C °^°°M (33) 

P{C DL ( A)= v 0i (A)} = f Cb(0 ) L 1 (A)-“ ,a »(l- ;PM (A)) C '’ ,0) -''“ (4) (3.4) 

where 


PLo(A) = -^(l-e-^)=l- Pu (A). 


(3.5) 
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Recall that in Ground Truth C Red Subunits are hypothesized to be in motion on the arc 
for time T; assume for now that T is long enough so that departure from the arc does not 
occur. These same transition probabilities prevail at every future PRT, i.e. at kA, k=l, 2, 
3,so long as the Subunit remains on the arc. 

It becomes necessary to specify Blue’s perception at t = 0 of the Red force size 
actually departing. We state some options for both observations and inferences. Not all 
are totally satisfactory. 

Option 1: Simple Moment Estimators 

It is easy to see from (3.1) - (3.4) that if C units transit 

£[c4a)|Q>(0)] = C D (0)p DD (A) + C L {0)p LD (A) 

(3.6) 

= Cd (0)pdd (A)+[C— Co(0)]/»i£)(A). 

Thus if C D ( A) is observed, giving Cd(A), a straightforward candidate moment estimate 
for C, is based on only the first two observations 

£(a) # = C °( A ) + Cd ( A ) ~ Pdd ( A )] . ( 3 ~ 

Pld{&) 

by stationarity an estimate based on any two consecutive observations is 

C(kA) # = + c -p((^ ~ 1 ) A )[/ ? ^( A )~ Pdd(A)\ 3 

Pld( A) 

Unfortunately these estimates, while technically unbiased, may assume negative values 
for certain observational values. To avoid this inadmissibility they may be modified to 


C(kA) = maxj^C(M) # ,o| 


which of course makes the estimate non-linear and hence awkward to study 
mathematically. 
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Time-Averaged Moment Estimates 

It is intuitively clear that, so long as Subunits remain on the arc, an improved estimate 
should result from an appropriate averaging process. One such is to simply average all 
past estimates with equal weights: if PRT is currently jA then 

C(jb) = -£c(kA) 

1 (3.10) 

= jc(;A)+^C((y-l)A) 

the latter recursive update formula is convenient. 

An attractive alternative is the exponentially weighted moving average (EWMA): 

C(yA) = aC{jA) + (1 - a)C((j -l)A), 0 < a < 1. (3.11) 

This estimate, e.g. with a= 1/3, automatically and mechanically reduces dependence 

upon the past in a way that the unweighted average, C, does not. It does so without 
regard to the relationship of the most current estimate, C(j'A ), to the immediate past, so it 

may be slow to adapt to events such as abrupt changes in arc occupancy, or to trends in 
time (adaptation comes by modifying a as a function of recent radical changes in C(jA); 

details remain for new work). But its simplicity and computational convenience cannot be 
denied or improved upon. 

Utilization of the Estimate in Model Context 

Within an evolution of the model it will be a simulated sequence of estimated 
C-values, i.e. either C(yA) or C(jA),j= 1,2,..., that will be presented to the analyst 

who guides one side’s actions. Corresponding estimates must be introduced into 
algorithms that govern the opponent’s behavior. These present estimates may be 
augmented by information as to the likely range of force-size (C-values) consistent with 
the estimates, and by further information concerning Subunit types; all such information 
elements or perception components are explicitly modeled as afflicted by (realistic) error 
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to some degree. The analyst will then make decisions concerning future moves, e.g. by 
Blue, in the face of the perceptions so delivered. 

Option 2: (Quasi) Likelihood Estimation 

Likelihood methodology is an alternative approach to combining information from 
previous observations. The awkward form of the likelihood function associated with the 
current Markov model for Co(kA) leads to use of a moment-matched Normal 
approximation or quasi likelihood, see Nelder and McCullagh (1983). Assume then that 
C(£A) has approximate conditional density 

expj - j k - (ck-\PDD + (C - Ck~\ )pld )1 (A j 

f{pk\Ck-\,C')~ --—--—--i- (3.12a) 

yj27t(J 2 k 

where 


- c k-\PDDi}~ Pdd) + {C- Ck-\)pLDi}~ Pld) (3.12b) 

and ck and Ck-\ are the observed counts at times kA and (k-l)A respectively. Now the 
likelihood function for unknown (but presumed constant) C is 

Z(Cj;data) = f[f{c k -,c k -uC ). (3.13) 

k =1 

Rearrangement puts this in the form 


j -^[ck-Ck-l(pDD-pLDyCp L Df / O 2 ^ 

L(C,j; data) = f]- r. ■ - 

*=i 


(3.14) 

-i[c 7 -Cj-X ( POD-PlD )-Cp LD ] 2 jo) j-\ -§C k -Ck-\(pDD-PLD)-Cpwf j 

= -- r _- Yl~ -r—_ 

Of course a\ depends upon the unknown C, but also upon all previous observations. We 
propose to iteratively update the approximate likelihood estimate of C as follows: 
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(a) write the likelihood up to (/-1)A in the normal form 


A[c-»i4hU 

L{CJ- 1; data) =- - . (3.15) 

^2xTj -1 


This replaces the product to j -1 on the right-hand side of (3.14); 

(b) complete the square to calculate the likelihood up to jA in the form 

L(CJ; data)= (3.16) 

■\j2,7r Tj 


where 


_ (mj-i / r y-i) + [( c y “ C J- 1Ota> “ Pld ))/ Pld]/(&]/Pld) „ ^ 

MJ - + l/(^L) (3 ' 17a) 

and 

r? = l/(l/r?-i + l/(Sj/plo)) (3.17b) 

where ctJ = c;.iPbb(1 - Pdd ) +max{(/rj_i -cj- 1 )p!.n(l - Pio),0} . 

Since g)-\ depends upon the unknown C-value, replace that value by 

Pj-\, available at (/—1)A; actually it is best to apply the formula (3.12b) altered to 
replace C - Cj.\ by max(/y_i-c 7 -_i, 0 ). This permits a numerical value for pj to be 
calculated; clearly this process may be iterated. 

(c) The above expresses uncertainty about C in a Bayesian style, effectively attributing 
to C a density with mean p.j and variance z] at time jA. Observe that pj actually is, 

as written in (3.17a), a linear combination or weighted average of the most recent 
moment estimator, C D {jA), as in (3.8), and the most recent point estimate, pj-\. 

The form is exactly analogous to the exponentially weighted moving average but 
with the weights determined in a systematic manner, i.e. updated by use of the 
estimated variance of the observations and the (approximate) model. Additionally 
there is the build-in assessment of the variability of the estimate furnished by r). 
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(d) The above setup can be made adaptive to changes, e.g. of arc content, by 
discounting the older likelihood. This leads to scaling up t)-\ by a factor K& > 1 : in 

(3.17) simply replace tj- X by Clearly more weight is now placed on the 

current estimate. 

Option 3: Maximum Quasi-Likelihood and Laplacianized Quasi-Likelihood 

An alternative to the simple iterative procedure of Option 2 is to collect all of the 
C-dependent terms of the quasi likelihood function (3.12) into the exponent: 

Q(c) = 

1 r / . , y,2 1 (3-18) 

--[ c k-{Ck-\PDD+\C-C k -x)p LD )\ - j -r—r--r---r 

2 Ck-\PDD\} - Pdd) + (C- c*_i )Pld\ 1 - Pld) 

- — ^\pk-\PDDi}- - Pdd) + {C-Ck- 1 )Pld{ 1 - Pld)\ 

Now minimize Q with respect to C so as to maximize the likelihood: differentiate with 
respect to C, set the derivative equal to zero and solve the resulting quadratic for the quasi 
mle; tedious algebra gives the result 


CMi.(k ) = 


1 

-b(k)+ jmax((b(£) 2 - 4a(&)c/(&)),o)j 2 
2a{k) 


(3.19) 


where 

a{k) = Pll(1~ PllY 


b(k) = \pll{\-Pll)\ +2(1- Pll) 2 cic-iPdd(1-Pdd) 

x(k ) = c k -c k -i[p DD -(l-P ll)] 

d{k) = -[2(1 - Pll)cic-iPdd{ 1~ PDD)*{k) 

+Pll (1 - Pll) x(kf - p LL (l - p LL )c k . x p DD ( 1 - p DD )]. 


The second derivative can be calculated and hence the approximate variance is 
<x 2 ml (k) = v(kf x [(1 - p LL f v(k) + 0 5[p LL (l -p LL )] 2 


(3.20) 
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where v(fc) = p LL ( 1 - Pu)C ML {k) J rc k ^p DD ( 1 -• 

If CMi{k) is now treated as normally distributed with the above parameters C ML (k) 
and (J 2 Mi(k) its contribution can be incorporated into an updating procedure analogous to 


that of (3.14)- (3.17). 

An alternative to the above procedure is to treat the quasi likelihood as a density for 
C, treated as a random variable a la Bayes. Rewrite (3.12) as 


(C-C)/(aC+b) 


g(C,Am)-K 


(3.21) 


where C is an abbreviation for Cd(A) # or Cd(&A) # as in (3.7), and a and b depend on 

data ck and q t_i, and parameters pdd and pld- 

(a) normalize, by determining K so that the approximate density 

« /( aC+b ) 

7 7 dC = 1 (3.22) 

J 0 4l7t4aC+b 

approximately, using Laplace’s method; cf. Tierney and Kadane (1986); do so by writing 
the integrand as eQ(Q, with 


2(C) = - i (C - C)*/(aC + b) - j lnfoC+ b) (3.23) 

and determine the minimizing value of C and the second derivative at that point; this 
leads to 


® -i( C-Cml) /o\ti 

K\ - - 7 =- dC=\ (3.24) 

J * r) 'rr 


from the previous development. 
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Similarly, 


E(aC+b ) = J (aC+b)g(C; data)^C (3.25) 

o 

can be approximated using Laplace’s method by writing the integrand as a Q m (Q with 


Q m (C) = -±(C-C m f/(aC+b)+^\n(aC+b) (3.26) 

and determine the minimizing value of C m , C m and the second derivative at that point. 
This leads to 


where 


E(aC+b) * ct 2 ml 

= L m 


<j 


2 

m 



In a similar vein 


(3.27) 


(3.28) 


E[(aC+ 6) 2 ] = J (aC+ bfg(C; data)</C (3.29) 

0 

which can be approximated by writing the integrand as e & '^ with 

Qy{C) = --(C- C v ) 2 /(aC+b)+-)n(aC+b) (3.30) 

2 2 

A 

and determine the minimizing value of Cy, Cy, and the second derivative at that point. 
This leads to 


where 


E[{aC+bf\ * Ke Qy{dv) ^ = 
= Ly 


(3.31) 


27 



cry =-1 


v_ 

a : 2 


- 1-1 

Qb'(G') 


(3.32) 


An estimate for C can be obtained from (3.27) 

A — b 

Cl =-. 

a 

An estimate for the variance of C can be obtained from (3.31) as follows. 
L v ~ z[(aC +bf ] = a 2 e[c 2 ] ■+ 2 abE[C]+b 2 . 

Thus, 


(3.33) 




L v - 2abCi - ft 2 


a 


Far[C] * ([iV -2abC L -b 2 ]/a 2 )-C\ 


(3.34) 


If Cl is now treated as normally distributed with the above parameters C £ and 
variance Var[C] its contribution can be incorporated into an updating procedure 
analogous to that of (3.17). 

Option 4: Numerical Integration of a Posterior for C Based on Quasi-Likelihood 

An alternative that is of strongly Bayesian flavor begins with the quasi-likelihood 
expressed as a density for C, i.e. (3.21), and numerically integrates to find E[C\ck, 1] 
and Var[C\ck, c*_i]. Following this, it performs the iterative update a la (3.17). These are 
the steps: 


(a) normalize (3.21), i.e. determine K so that 


4 


q ^2^(ax + b) 
where \(ax + b)=\ iiax + b>Q and is 0otherwise; 


\(ax+b)dx = 1; 


(3.35) 
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(b) numerically integrate to obtain 


;) /(ax+b) 

-- — 1 [ax+b)dx; (3.36) 

ax+b) 

and 


(c) numerically integrate to obtain 

«> ~{x~cf j{ax+b) 

E\c 2 \c k ,c k -A = K\x 2 e ■ .■■■■■ ==- 1 (ax + b)dx. (3.37) 

0 yj2x(ax+b) 


(d) Convert the latter numbers into a value for the Var[C]. 


(e) Apply (3.17) to update. 



4. Sensor Models Revisited: Binomial, Multinomial, and Generalized 

This section reports models for describing the output of imaging sensors in an 
operational context. They are applied in certain simulation (Monte Carlo) models 
(JSTOCHWARS/JWAEPS) to represent detection and (mis) classification of entities 
(enemy vehicles) on the surface of the earth from points above it, where these observation 
points can either be fixed in place (hovering) or in motion. 

Basic Idea 

A sensor’s view of the earth surface is limited to a spatial region thereon; call this the 
footprint-, in JS this is over an arc, or a node, but in reality is a (two-dimensional, 
sometimes more) region. The footprint can be virtually stationary, or can move in search 
of new opportunities. An image of the footprint elements detected by the sensor is 
portrayed on a screen. It is initially assumed here that, if there are C (generic) items 
(potential targets) on the surface (such as tanks, other wheeled vehicles, trucks) detectable 
by the sensor system, i.e. in motion for some (e.g. J-STARS), then these are seen on a 
glimpse or scan with conditional probability p (and not seen with probability q- 1 -p); 
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the events of seeing on successive scans are conditionally independent. Alternatively, a 
picture is on the screen, but the scan or glimpse is a time of residence of an operator 
foveal image over a subregion of the screen. Items within that subregion are candidates 
for operator attention (mis)identification, and communication towards weaponeers. The 
physical conditions that affect the probability of detection are (a) terrain, or sea state, 
including vegetative coverage, (b) the (radar) cross-section of vehicle types of interest; 
same for infrared, heat detection, (c) state of motion, or hiding, (d) action by enemy (Red) 
entities, such as combat/attack by them that facilitate detection, such as SAM site 
radiation during air defense (AD) activities. 

Models 

An initial model is that R, the random number of entities (e.g. potential targets) 
revealed or detected on one sensor scan/glimpse is binomially distributed, conditionally 
on environmental factors. Likelihood and moment estimators approximately provide this 
simple estimator of C, 



(4.1) 


where • signifies explanatoiy variables. A candidate parametric model for this 
dependency is the random hazard (or variations thereof, to be described); see Gaver 
(1963), 


p(v,s) = e-<- hK ^ (4.2) 

/ 

where In /? 0 (v) = ^ , a regression term. 

i=0 

The vector variable v = (vi, V 2 ,..., v p ) is one of observable explanatory variables such as 
range from sensor to target region, terrain type, cross-section of targets, atmospheric 
properties, etc. The (possibly vector) variable e= (e u e 2 ,..., e m ) describes random 
environmental fluctuations; these may be common to several, e.g. consecutive, sensor 
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scans because of cloud cover. The coefficient vectors J3 can be estimated from data, and 
modified for sensitivity checking purposes. The idea is that a sensor output is modeled as 
the outcome of C biased (not necessarily !4 - 14) but independent coin flips with 
probability of success simultaneously influenced by specified explanatory variables. The 
advantage of such a regression model is that it can be fitted using observations made 
under different conditions, i.e. on different target types, with different ranges, terrain 
types, i.e. by pooling data to estimate the regression coefficients J3. There are standard 
computer programs for carrying out maximum likelihood calculations. Then the 
expression, e.g. (4.3) can be used to predict the actual number of targets present. 

Notice that the estimate now has the form 


C(v,s) 


R 

p{v,e) 


(4.3) 


This is conditional on e, so this model has the effect of increasing the observed number, 
R, to account for those targets not observed; of course this adjustment is model- 
dependent. 

Since, conditional on e. 


E[R | s] = Cp(v, s), 
the unconditional expectation of R is 


(4.4) 


E{E[R\e]} = CE[p(v,e)]; (4.5) 

so, if we assume the distribution of p(v, e), given the explanatory variable values, y, is 
well-enough known, an unbiased (approximately!) estimate of C is 


E[p{y,s)] ' 

To obtain a numerical estimate, having observed R at t, R t = r t , we write 


(4.6) 
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n 


(4.7) 


E[p{v n e t )] 

where it is being assumed for the present that { s t } is a sequence of independent random 
variables with “known” distribution. 

Note that we are allowing here (modeling) a situation in which all observers are 
simultaneously affected by a chance event, such as cloud cover. We are not (yet) 
representing adaptive action by the sensor-carrying platform, such as flying below the 
clouds (not always possible), or using different sensor technology, e.g. IR. 

Variance 

The estimate at time t, given by (4.7), is an instance of a random variable, (4.6). 
Conditional on e t , and using the binomial model, 



_ Cp[v t ,s t ') 
E[p{v t ,s t )} £[/?(y,,£,)] 


(4.8) 


whose expectation is the true value, C. 

Consider the variance of C t , a measure of the variability of C, around C, the true 
value. .Note that this assumes just one glimpse “at C” while a given e, is present, if e t 
remains the same for several glimpses the model must reflect this, and can. First, 



C/?(y,,<g,)(l-.p(y „£t)) 

(£[/?(y„f,)]) 2 


It is known that, unconditionally, 

Var[c t ] = E{Var[c t |*,]} + Var e , {E[c t \s t ]} 

_ c [ E [p{y.t ,*;)]- e [p 2 (k, > £t )]) + C 2 Var[p(v t , s t )] ( 4 - 9 ) 

[E[p{v t ,s t )]f 

In order to compute an estimate of the variance replace C by C as in (4.6) or (4.7). 
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Analytical Model for Detection Probability, p(v, s) 

A candidate, and convenient, parametric analytical model is the following; described 
in (4.2): 

p{v,s) = exp[-^(-lnp 0 (y))] (4.10) 

where we can let 


-lnp 0 (y) = exp 




= exp[/? 0 +/?E • 


(4.11) 


L i=l J 

Note that one can explicitly evaluate all expectations in (4.9) for distributions whose 
Laplace transforms are available; this is a wide class that includes the gamma, positive 
stable, inverse Gaussian, and many others, including convex mixtures of the above. We 
give only the gamma example. 

Gamma. Here 



and£[f] = 1, Var[e] = 1/a. 

The correction factor for the mean in (4.7) is seen to be 



(4.12) 


(4.13) 


which, for a= 1, is just a logistic regression function; other values of a>0 simply 
generalize the model, possibly usefully. By simple analogy terms occurring in the 
variance are evaluated 


4> 2 M=y 


1+ ! e (p° + i *) 2 

a 


(4.14) 


Var[p{s)\ = £[p 2 (y,f)]-(£[p(y,^)]) 2 can be evaluated using (4.14) and (4.13). 
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Explaining Variability Alternatively 

The previous model allowed for variability in detection probability to be represented 
by deterministic explanatory/regression variables, v, and the initial binomial model to be 
extended by additionally randomizing simultaneously all conditional detection 
probabilities on a scan/glimpse. This device can be used to represent the effects of 
simultaneous random environmental variations unexplained by regression. 

An alternative possibility is to explain variations between potential targets (e.g. 
because of different orientations or partial cover) by randomly introducing an s-effect for 
each target', an independent and identically distributed assumption is convenient and 
parsimonious. This approach will be pursued in later work. 

Now describe the occupancy of a node or arc, denoted by n, by the numbers of units 
of different types (e.g., heavy armor brigades, light armor brigades, etc.); {/,(«; t ) is the 
number of such units of type i (i = 0, 1, 2, ..., I) at time t at node n; Uo(n; t ) can designate 
the empty node if necessary. 

Within a unit of type i, there may be several asset types, such as tanks, armored 
personnel carriers, etc. Assets may also refer to subunits (sunits) such as tank companies. 

Agree that distinguishable assets occur in J classes, and that a unit of type i has a mean 
number of assets of type j (j= 1, 2, ..., J) equal to cijj(n,t), and a variance cr?j(a;n,t). The 

initial values of the mean number of assets of type j for a unit of type i may be taken from 
the Table of Organization and Equipment (TOE) for that type of unit. Furthermore, adopt 
the provisional model that the actual number of assets of type j possessed by a particular 
randomly selected unit is a random variable, A(i,j, n, t ) with distribution function 
Fij(x; n, t). When convenient, and for illustration, we take A$J; n,t ), k= 1,2,..., Uu 

i.e., the numbers of assets of type j owned by the 17, copies on Node n to be independent 
and normally/Gaussian distributed. The time-dependent parameters ajj and cr| can reflect 
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the fact that a campaign has been in progress for some time and attrition has occurred and 
is subject to change. Let 

_ u i 

A (i, j, n,t) = £ A ( (i, j,n,t) 

t=\ 

denote the total number of assets of type j possessed by units of type i at node n at time t. 

It is convenient to refer to the vector of distributions of typical asset-type numbers for 

a particular unit type as the signature (or asset signature) of the unit type. Note that 

signatures of different individual units of the same type will inevitably differ if their 

various asset counts differ, as could well happen. Let 

_ 1 u < 

A {j,n,t) = Yj X Ai *) 

i=1 1=1 

the total number of assets of type j at the node at time t. 

Finally, let Sj(s, n; t) denote the total count of assets of type j by sensors of type 5 
(s = 1 , 2,..., S) at node n at time t. Sj represents the quantitative perception of the 
opponent’s type j asset level. 

In this section and the next two, we present and more thoroughly investigate models 
for sensor observations of assets present on nodes or arcs that are based on the binomial 
and multinomial distributions. We consider the behavior of approximate procedures for 
updating the estimate of the number of units or assets on the node/arc using the binomial 
and multinomial observation models. The procedures are based on approximating the 
sampling distribution of the estimate of the number of units by a normal distribution. The 
paper by Hall (1994) gives some insight into when these approximations may be 
appropriate, but their convenience for the purpose of minimizing computation time is of 
overriding concern in the context of a computer model. 
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4.1 Combining Binomial Observations 

A model for two independent observations of a unit with r assets by two different 
sensors is as follows. Let X\ and Xj be independent random variables having binomial 
distributions with common number of trials r and success probabilities p\ and pi- The 
problem is to estimate r from observations x\ and xj. 


4.2 The Likelihood Approach 

The likelihood function is 


L(r,x u x 2 ) = Y\ 


2 (r ^ 


i=i 


W 




An approximate maximum likelihood estimate for r can be obtained by solving 


T(r;x i,s 2 ) =1 
L(r-V,x x ,x 2 ) 


for r. This results in the equation 


(4.15) 


(4.16) 


(1-a)—— (1-^ 2 ) = 1- (4.17) 

r-X\ r-x 2 

Thus, r satisfies the quadratic equation 

r 2 [l-{l-p l ){l-p 2 )]-r(x l +x 2 )+x ] x 2 =0. (4.18) 

If there are 3 observations a similar argument would result in r satisfying a cubic 
equation, etc. We get an (approximate) maximum likelihood estimate, with the good 
properties of the m.l.e., but must solve an awkward equation. This may not be a problem 
for a small number (especially one or two), but if many must be done quickly, it may be. 
We look for an alternative. 


4.3 A Moment Estimator and the Normal Approximation 

The moment estimator for r is seen to be 


36 



(4.19) 


Note that 



rf-l r iPi 
E[ri\ = — = n 

Pi 

and so n is unbiased, as it was designed to be. Further, 

Var[r,] = -^Va^X,} = IBEzHI = HUzIil . 

Pi Pi Pi 


(4.20) 


(4.21) 


Jt* 

Approximate the distribution of tj by a normal distribution with mean ju(i) = — and 

Pi 

variance v(z) 2 = x t (l -pi)/pf ; if x z = 0 then set v(z') 2 = vo > 0 where vo is chosen by the 
analyst. To obtain a combined estimate of r from f\ and r 2 take the weighted average 
with weights the inverse variances: 


The estimated variance is 


n *2 

v(l 

1 1 

v(1) 2+ v(2) 2 


(4.22) 


a 2 = 


1 1 


v(l) v( 2 ) 2 


(4.23) 


Figures 15-17 present results of a simulation experiment with 499 replications. For 
each replication two independent binomial random numbers are drawn, each with 20 
trials and one having probability of success 0.5 and the other having probability 0.9. In 
each replication, the estimate obtained by solving (4.18) is computed and the estimate 
(4.22) is computed. 

Figure 15 displays a histogram of the differences of the likelihood estimate and 
moment estimate. Note that the likelihood estimate tends to be somewhat larger than the 
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moment estimate. Figure 16 (respectively 17) displays the histogram of the likelihood 
(respectively moment) estimates. The moment estimator appears to give satisfying results, 
although likelihood is probably superior. 


4.4 Approximate Normal Updating Using Binomial Observations 

Let the current estimate of the number of subunits or assets on a node, r, have a 
normal distribution with mean p t and variance a ]. Suppose a new binomial observation, 
Xf+\, arrives with parameters r trials (there are, and have been, r units present) and 
probability of detection p. We assume the updated estimate has a normal distribution with 
mean 


Mt , x t+ i/p 

hr + i=-^-(4.24) 

of v,+i 

and variance 


where 


_2 _ 
07+1 - 




(4.25) 


V/+1 


_J*m(1 ~p)Ip 2 


Vo 


if x t+ \ > 0 
if*/+i =0 


(4.26) 


where v© is a constant chosen by the analyst. 

Note: the above algorithm is only appropriate if the latest observation is of (about) the 
same number of units as the earlier ones: not much has changed. But changes will occur 
and must be detected so the above is only a beginning. 
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5. Beta-Binomial Observations and Their Combination 

One generalization to the simple binomial model of Section 4 is to allow the 
probability of detection p to be random with beta distribution having density function 


Ap) = 


P a -\1~P) 


P -1 


B{a,0) 

0 


if 0 < /? < 1 
ifp < 0 or p > 1 


(5.1) 


where 


B(a 


The beta random variable has mean m = 


a 


a + P 


(5.2) 

and variance 


vl-apj (a+0) 2 \(cc +/?) +1] . The purpose of this generalization is to recognize 


variation in the actual probability of detection,/?: make p a random variable. 

Let X be a random variable whose conditional distribution given p is binomial with 
parameters r trials and probability of success p; now additionally suppose p has a beta 
distribution with parameters a and p 




r(a + 0) r(a + x)r(N - X + p) 
r(a)r(/3)r(N +a+ /i) 


(5.3) 


An estimator for r is 


r = 


X 


m 


E[r] = = ^E[rp] = r. 


(5.4) 

(5.5) 
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Thus, r is unbiased. Next, evaluate the variance: 


Var\f~\ = + Far^£[r|/>]J 



rp(\-p) 


+ Var\ 



_r _ afi _ r 2 af3 

m 2 (a + f} + \)(a + 0) m 2 (a+/?) 2 (a + /?+l) 


r aft 

m 2 (a + 0){a + /?+ 1 ) 


1 + 


r 

(a + P) 


In summary, 


(5.6) 


Var[f] = ~vl[(a + p) + r] = + P + r}. 

Notice that this estimate depends quadratically upon the unknown r, rather than linearly 
as in the plain-vanilla binomial case. 


5.1 Approximate Normal Updating Using Beta-Binomial Observations 

Assume the current estimate of the number of subunits on a node, r, has a normal 
distribution with mean and variance cr 2 . Suppose a new beta-binomial observation, 
Xf+\, occurs. 

The estimator for r using X t +\ is r t +\ = X t+ \/m and we approximate its variance by 
Vf+i = f l+ 1 fi[a + f3+r,+ 1 ]/ a(a + fi+1 ) if x t+ \ > 0 and v 2 +1 = vo if*/+i = 0. 

We assume the updated estimate has a normal distribution with mean 


Mt+\ = 


, r l+ \ 

2 ' 2 

q_i v<+i 

1 1 

2 + 2 

<?t V,+i 


(5.7) 


and variance 
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__2 _ 
&t +1 - 


1 


1 


1-1 


2 2 
Vt V /+ i 


(5.8) 


5.2 Maximum Likelihood Estimate for the Number of Trials, r, in which the 
Observations are iid Beta-binomial 

Let X\,..X t be independent beta-binomial random variables with parameters r trials 
and beta parameters a and (3. The likelihood function for r is 


L I (r;x u ...,x t ) = K— 


(rif 


fjr(r-x/+yff) 


1=1 


n(r-*,) ! r ( r + a + / 3 )‘ 


(5.9) 


1=1 

The maximum likelihood estimate is that integer r which maximizes Lj. An approximate 
solution can be found by determining that r such that 


Lijryci . Xt) 

L I (r-l,x u ...,x t ) 


For t ~ 2, for r > max(xi, X 2 ) + 1 


J, r (/-;jci.jc f ) _ r 2 (r-xi+fi-j)(r-X2+fi-l) ^ ^ 

L 1 {r-l’,x\,...,x t ) {r-x l )(r-x 2 ) ( r-l + a + /3) 2 

Setting the ratio equal to 1 results in a cubic equation for r; the minimum solution larger 
than max(Ai) is selected as the reasonable point estimate. 

Figures 18 and 19 display results of a simulation experiment with 99 replications. 
Each replication consists of 2 independent beta-binomial random numbers; the binomial 
number of trials is 20 and the parameters of the beta distribution are a= 5 and (3=2. 

Figure 18 (respectively Figure 19) displays a histogram of the moment estimates 
(respectively the maximum likelihood estimates). Note that the moment estimate is much 
more variable than the maximum likelihood estimate. 
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5.3 Common but Random Probability of Detection 

Assume X\, ...,X t are conditionally independent binomial random variables given p 
with r trials and probability of success p but p has a beta distribution with parameters a 
and /?. Randomization of p common to all sensors could represent a common effect on 
them all, e.g. by weather, or cloud cover, or perhaps terrain variation. Then 

P{X i =x u X 2 =x 2 ,...,X,=x t } 

a+ S^l r f^ + E( r -^) 

1=1 l V *‘=i _ 

;.!UJJr(a)r(/?) r(«+/?+»•) 



(5.11) 


Thus, the likelihood function for r is of the form (K an arbitrary constant) 


, f r N r ^+Z( r -^) 

L{r,x x ,...,x t ) = KY[ -t- 

T(a + p+tr) 


(5.12) 


For r > max(x,) +1 and t = 2 


Ljr^x 

L(r-l;xi,...,x t ) 

(5.13) 

r 2 (r-x\+r-x 2 + /3-l)(r-X\+r-x 2 + fi-2) 

(r-x\)(r-x 2 ) (2 r + a + fl-l)(2r + a + fi-2) 

Setting the ratio equal to 1 results in a quartic equation for r which can be solved 
(numerically); the minimum solution greater than max(jc,) is the reasonable one. 
Numerical stability concerns suggest that the logarithm of the likelihood function be 
computed and Stirling’s formula (cf. Feller [1968]) be used to approximate the gamma 
functions; that is, compute 
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' (r\ r ii r i 

^(r;*,,...,*,) = £ln + p+^{r- Xi )-j- In/?+^(r-x ; )-l 

;=i \ x jJ Li 2 J L 

~ P+Y,(r-Xi)-\ 

i 

a + fi+tr-^ ln[a + fi+tr\ + [a+fi+tr-\\ 

= E hi f ^] + P + y L ( r ~ x ‘)-\ 111 P + J j ( r ~ x i )- 1 

j=\ \ x jJ Li z \ L / 

1 

- a + B+tr — 

L 2 

where K does not depend on r. 

Figures 20 - 23 report the results of simulation experiments. For each replication one 
random number, p, is drawn from a beta distribution with parameters a and ft. For the 
experiments reported in Figures 20-21 (respectively Figures 22 - 23) two (respectively 
10) independent random numbers are then drawn from a binomial distribution with 20 
trials and common probability of success p. Conclusion: for cases considered, the 
moment estimator appears to be an adequate approximation to the maximum likelihood 
estimate. Note that, as mentioned earlier, the above Beta binomial model represents extra 
variability common to sets of binomial probabilities (of success). The present Beta 
binomial does not represent any extra binomial randomness between the outcomes of 
successive looks at the same unit by the same sensor (within variability). 

6. The Many-Unit Type (Multinomial) M^classification Problem 

We next present a more realistic version of the models considered earlier in this 
report. In an actual operating environment there is presumed to be a melange of different 
observables: members of several classes of different subunits, such as tank companies: r z - 
subunits of type i, i= 1,2, ...I in a specified part of the theater (on a given arc or at a 


ln[a + (3 + tr] + K 
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given node). Assume that each subunit type has a probability of being detected on a 
sensor pass or observational opportunity; let this be dj for subunits of type i. Furthermore, 
once detected the subunit is assumed correctly classified with probability ca, but 
incorrectly so, and to one of the recognized classes, with probability Cy, j * i. It is 
optimistic, but convenient, to assume that cu « 1 and cy, j * i, is close to zero; this will 
not always be necessary but can be a useful start at times. An overly simplistic way of 
handling the misclassification problem is to simply assume that the probability of 
misclassifying i to j * i is a small constant. 

Note that each of the above parameters is likely to be specific to sensor type: di(s), 
letting the sensors be classified as of type s = 1,2,..., S. And, during a particular 
(e.g. 2-hour) reference time period the number of opportunities that a given sensor has for 
“seeing” a unit of type i at a given location will vary depending upon the action(s) of the 
subunit(s) and the number of passes the sensor makes over the particular region under 
examination. This will be a programmed quantity and will be known. The actual 
detectability of a subunit, i.e. as influenced by the subunit present where looking occurs is 
unknown and must be estimated. 

6.1 Model 

Begin by considering one pass by a sensor over a region (e.g. arc). First, carry out 
detections: let D z be a random variable denoting the number of detections of (sub)units of 
type i; if desired specify Dj(s). Assuming sampling with replacement is adequate, 

Di is Binomial (r z , dj). (6.1) 

We assume that Dj is not directly observable. Of course there is identity information 
inherent in what different sensor types see: if sensor type 1 can only see subunit type 1, 
sensor type 2 can only see subunit type 2, then if d\ - 10 and di = 0 it is clear what types 
of subunits are present. 
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Next, classify all of the detections according to the probabilities cy. Consequently 
each of the D,- (anonymous) observations is independently identified with some class j, 
where j= 1,2,Thus r/, the unknown number of type i subunits gives rise to a 
multinomial distribution of numbers of values Xy 

P {•X’ii -xt\,X l2 =Xj2,...,Xu = xu } = p*f (l- d t ) r ' (6.2) 

■ Xu J j=i 

where py = djCy. 

But what can be observed is, say, 

/ 

x j =Y J x i j (6.3) 

i=i 

all observations classified as type j. The objective is to turn the above into a reasonable 
estimator of 77 , assuming that di and cy are known, at least initially. As stated before, 
there is useful classification information in xj(s), 5 = 1,2,..., S if sensors are differentially 
sensitive. 

Moment Estimator 

Apparently a relatively quick way of estimating the parameters r\, r 2 , is to use 
the method of moments. This may not be as efficient as other ways, but conveniently 
makes low computational demands. 

Let Xy{s) denote the number of times a sensor of type 5 sees a unit of type i and 
classifies and reports it as of type j. Then 

E[Xij ( 5 )] = rfit (s)cij ( 5 ) ■ npy ( 5 ). (6.4) 

/ 

Now Xj(s) = £x,j(s) models the actual observed data xj(s), so the method of moments 
1=1 

puts 

/ 

*/( 5 ) = &#/( 5 ) 7 = 1,2,...,/. (6.5) 

1=1 
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In matrix terms 


so, if the inverse exists 


x(s) = rpO) 


( 6 . 6 ) 


r = x(s)p _1 (*) (6.7) 

This is directly analogous to the previous scalar setup. But a solution of / simultaneous 

/ 

equations is now required. Notice that ^Pij(s) = d t (s ) < 1. 

j=i 


Approximations 

The solution offered above may be statistically somewhat inefficient (use of moments, 
rather than likelihood or Bayes) and is also computationally troublesome. Here are some 
approximations; these are tentative and will be checked by simulation. They are based on 
the assumption that ptfs), j * i, can often be expected to be much smaller than pu(s'). If 
not, then we are sometimes led to say, skeptically, “I have just seen a goat, but I know 
that the other side always dresses sheep like goats, so I will treat this as a sheep.” Such is 
related to decoy and deception issues, but also to imperfect Battle Damage Assessment 
(BDA), and to occurrence of false targets in general. 

Assuming the above to be true we can plausibly estimate r, initially by the first 
(naive) estimate : 


Pii\ S ) 


( 6 . 8 ) 


The model for this is 



Now take expectations: 



Now since p/d «pa the second term can well be small (unless r/-» r,). So a possible 
way to correct (and improve) the initial estimator is to approximately remove the 
estimated bias, i.e. compute the second estimate'. 


rP(s) = max . (6.11) 

V k*i Pii\ S ) 

It might be worth iterating this to obtain rp\s), etc., but this step will not be taken. To 
obtain an estimate for the variance of f£ 2 \s ), note that 


Var[Xj ( 5 )] = £ r iPij (s)[l - Pij ( 5 )] 

1=1 

which can be estimated as 


Var[Xj ( 5 )] s Yj 9 i 2) 0 s)pij (*)[l " Pft ( J )] 

1=1 

for example. 

Approximately, 



( 6 . 12 ) 


(6.13) 
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if A/(5) - 0 set Var[Xj(s)] = vo where vo is chosen by the analyst. Of course this last is ad 
hoc, its properties are best understood by use of some simulation. 


6.2 Approximate Normal Updating Using Multinomial Observations 

Let the current estimate of the number of subunits of type /= 1, ...,/on a node, 77, 
have independent normal distributions with mean /u t {i) and variance a 2 t (i ). Suppose a 

new multinomial observation X t +\ arrives from a situation with parameters 77, pij(s), for 
i= l,I and 7 = 1 As usual 77 is unknown, but py(s) is assumed known. 

We assume the updated estimators have independent normal distributions with means 




MO , h 

< 7 j{f) vf+i(i) 
1 1 
cr 2 t{i) Vf+i (7) 


and variance 


(6.14) 


of*i(0 = 

where r t is an estimate of 77 and vf+\ (z) is 


—-p" (6-15) 

of(0 v, 2 +1 (z) 

an estimate of the variance of r, . The estimate 


of the variance uses ( 6 . 12 ) - ( 6 . 13 ). 

A Simulation Experiment 

In this section results of a simulation experiment are reported. 

In each replication of the simulation two multinomial random numbers are generated. 
Each has number of trials (r\, 7-2), probabilities of detection (d\, c/2) and probabilities of 
classification ( c\\, c\2, 021, c 22 )- The observed data are 


A 1 (l;s) = A n (l^) + A 21 (l;5);A 1 (2; J ) = A n (2;5) + A 2 ,(2;5) 


and 
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X 2 {l;s) = Xn{l;s)+X 22 {l;sy,X 2 (2;s) = X l2 (2;s) + X 22 {2;s ). 

Two procedures for combining the two observations are considered. In one, called 
Average, the average of the two observations is computed and the estimate (6.11) is 
computed using the average X\,X 2 where 

*i=[*i(l ;s) + Xi(2;s)]/2. 

Figures 24 and 25 present histograms of the moment estimate using one observation 
(Xi(l; s ), X 2 (l; s )) obtained by solving simultaneous equations (6.6) and histograms of 
the approximate moment estimate (6.11). The two estimates appear to be reasonably close 
with the approximate moment estimate being slightly lower. Table 6.1 presents statistics 
for the simulation results. 

Table 6.1 
1 observation 

d\ = 0.7, di = 0.9, cjx = 0.8, ci 2 = 0.2, c%\ = 0.1, C22 = 0.9 
500 replications: ri = 10, r-i = 15 


Proc 

Mean 

of Estimates 

Variance of 
Estimates 


A 

n 

A 

ri 

-A A 

n r 2 

Simul Eqtn 

10.3 

14.8 

14.8 7.34 

Approx 

10.0 

14.4 

14.0 6.9 


Figures 26 - 27 present histograms of the estimates resulting from the approximate 
normal procedure for combining two multinomial observations (Xi(l; s), ^ 2 ( 1 ; 5 )), 
(X\(2; s), XjQ.', 5)); the estimates of r[ combined are those resulting from (6.11), (6.12), 
(6.13), (6.14), (6.15). The figures also present histograms of first averaging the two 
multinomial observations and then applying estimator (6.11)- (6.13). The two 
procedures appear to yield comparable histograms. The estimated variances are computed 
using (6.13). The mean of the estimated variance is always smaller than the variance of 
the estimates. Table 6.2 presents statistics of the simulation results. 


49 



Table 6.2 
2 observations 

d\ — 0.7, di = 0.9, c\i = 0.8, cj2 = 0.2, cji = 0.1, C 22 = 0.9 
500 replications: ri = 10, rj = 15 


Proc 

Mean 

of Estimates 

Variance of 
Estimates 

Mean of 
Variance 
Estimates 


h 

h 

h 

h 

f\ h 

Ave 

9.8 

14.7 

6.1 

3.6 

5.8 2.8 

Normal 

9.4 

14.9 

6.7 

6.0 

5.7 4.9 


Table 6.3 records statistics from a simulation with 500 replications. For each estimate in 
a replication the following confidence interval is computed: ^ - 2-^Var , r, + 2^Varr t j; 
this confidence interval should have coverage at least 95%. Table 6.3 records the fraction 
of confidence intervals which cover the true value of Note that the fraction tends to be 
a little smaller than 0.95. 

Table 6.4 reports results of a similar simulation but one in which r\ = 5 and ri = 8. 
Note the fraction of confidence intervals that cover the true rj tend to be smaller than 
0.95. This behavior is due to the greater likelihood that the estimate of is 0. 
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Table 6.3 

Statistics of Simulation 

d\ = 0.7, di — 0.9, c\\ = 0.8, C 22 = 0.9, c\2 ~ 0-2, C 21 = 0.1 
500 replications 
r\ = 10, 7*2 = 15 


No. 

of 

Obs 

Type 

of 

Unit 

Estimator 

Mean of Est 
(Std Dev) 

Mean of F 

Est Std Dev 

faction of inti 
n ± 2 *JVar Pi 

covers true v 

srvals 

that 

alue 

1 

1 

Approximate 

9.8 


0.93 



Moment 

(3.5) 

(3.4) 



2 


14.8 


1.0 




(2.7) 

(6-6) 


2 

1 

Normal Updating 

9.6 


0.93 



with Approx 

(2.5) 

(2.4) 




Moment est. 





2 


14.4 


0.94 




(2-3) 

(2-2) 


2 

1 

Average obs. first 

9.9 


0.96 



then approx 

(2.4) 

(2.4) 




Moment Est 





2 


14.5 


0.91 




(1.9) 

(1.6) 
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Table 6.4 

Simulation Statistics 

r\ = 5, ri = 8, d\ = 0.7, di = 0.9, cn = 0*8, £22“ 0.9 
500 replications 


No. 

of 

Obs 

Type 

of 

Unit 

Estimator 

Mean of Est 
(Std Dev of 
Est) 

Mean of 
Est Std 
Dev 

Fraction of 

- 2 y]Var r ( , %+ 2^jVar r, j 
that covers true r,- 

1 

1 

Approximate 

5.0 


0.93 



Moment 

(2.6) 

(2.5) 



2 


7.8 


1 




(1.9) 

(4.7) 


2 

1 

Normal Updating of 

4.6 


0.87 



Approx Moment est. 

(1.9) 

(1.7) 



2 


7.8 


0.94 




(1.7) 

(1.6) 


2 

1 

First average obs. 

4.9 


0.91 



then use Approx 

(1.9) 

(1.7) 




Mom. Est. 





2 


7.8 


0.91 




_ US _ 

(1.2) 



Observe that the coverage fraction in the above table is nearly always too low, i.e. is not 


conservative. Practically, one might simply change from a multiplier by 2 to something a 
bit larger, e.g. from the Student’s t distribution. 

All of the above could be redone in a Beta-mixed context. 

7. Some Approaches to Updating Perception At a Node Using 
Information from Neighboring Nodes 


7.1 Introduction 

Describe the occupancy of a node or arc by the number of units of different types that 
occupy the node or arc; the node or arc may well be empty. Since units move along arcs 
and nodes, the occupancy of a node/arc changes over time. Each unit may have several 
types of assets, e.g. tanks, soldiers, APC’s (Armored Personnel Carriers), etc. and 
different types of subunits, e.g. tank companies. A protagonist gains information about 
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the node occupancies through the use of sensors which give him error-prone observations. 
In an idealistic and theoretical sense it would be desirable for a protagonist in a theater- 
level campaign to possess the joint distribution of the types and numbers of units 
occupying all relevant nodes, along with their assets, at every node for all times t. This 
distribution will be influenced by the protagonist’s allocation of sensor assets, by his own 
force maneuver as well as by the opponent’s maneuvers to accomplish a planned course 
of action. 

However, the development of an automated systemic model of a theater campaign 
must realistically stop short of a completely integrated joint probabilistic assessment of 
node occupancy; computational burdens are simply too great. 

Algorithms for updating perception of the numbers of assets/subunits and numbers 
and types of units at a node, using local information for that node, have been discussed 
previously; see Appendix A and Sections 4- 6. In this section we investigate the 
behavior of alternative perception updating algorithms. Some of these algorithms use 
information from neighboring nodes or arcs. 


7.2 The Model 

For simplicity of discussion and notation we will assume for the present that there is 
one type of unit, e.g. brigade, with only one asset type. Let K(n, t) be the number of units 
that occupy node n at time t and let A(n, t) be the number of assets or subunits at node n at 
time t. Assume the prior distribution for the number of assets given the number of units 
with no sensor observations is 


P{4«,0) <=da\K(nfi) -k}- expj ^ ( j P-l) 

where /is the mean number of assets associated with one unit and cr 2 a is the variance of 
the number of assets associated with one unit; y is derivable from the Table of 
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Organization and Equipment (TOE) for the unit plus any information available 
concerning unit campaign experience. 


We assume one kind of sensor which counts assets with 


P\S edx\A{n)-= a\ = -^=- exp 
1 J -y/2 TtT 


fx-a\ 2 


\ r 


(7.2) 


For illustration purposes suppose there are three nodes. 



Suppose it is known that units at node 1 will, after a random time T\, move to node 2 
and after a further random time 72, move to node 3. 

The problem is to use sensor observations to estimate the number of units and number 
of assets at each node. In general, there will not be enough sensors to observe all three 
nodes at all times. 

Let the marginal distribution of the number of units on node n at t, given all sensor 
observations up to t be 


7t(k;n,t ) = P^K{n,t) = Ar|all sensor observations during [0,t]j (7.3) 

and let the marginal distribution of the number of assets at node n, given sensor 
observations and units be 


dF k (a;n,t ) = P^A{n,t) eda\d\\ sensor observations during [0,f], K(n,t) = &} 

expj- ^ (a - m k (n,t)f /vf («,*)}• 


(7.4) 


42ftv k (n,t) 


In the remainder of the paper we present algorithms for updating n and Fk as more 
information becomes available. 
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7.3 Updating Using Same Node Information 

Suppose there is a sensor observation x(n, t + 1) of the assets at node n at time t+ 1. 
One can update the perception at node n using only the sensor observation at node n as 
follows. 

(7.5) 


v* (»,/ + !) = 


v k (n,t)/a(a) 


m k (n,t + l) 


m k {n,t) x(n,t + 1) 
v| (n,t)/a(a ) + t 2 
1 

v k (n,t)/a(a) + t 2 


(7.6) 


7r(k,n;t + 1) = Cn{k,n,t)- 


427t.lt 2 + 


= Cn{k,n,t)w[k',n,t,x ) 


AM 

a{a) 


exp* 


r 

- -i' 1 

1 

(x(n,t + l)-m k (n,t)f 

i CN 

T 2 + V|M 


a {a) j 


(7.7) 


where 0 < a(a) < 1 is a constant discount value. If a(a) = l, then the procedure is the 
same as that in Appendix A and Sections 4-6. 

If there is no sensor observation of node n at time t + 1, then set 


m k (n,t + l) = m k (n,t ) 

v| (n, t + i) = v 2 (n,t)/a 0 (a ) 
for a discount constant 0 < ao(a) < 1. 


(7.8) 

(7.9) 


7.4 Updating Using Neighboring Node Perception 

Since it is assumed from the combat situation of the section that the units that are on 
node 1 will eventually move to node 2, it would be advantageous to incorporate this 
information in the updating. 
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a. Modifying the Prior of the Distribution of Numbers of Units 

One consequence of a Bayesian procedure is that as more evidence accumulates for a 
particular number of units at a node, the procedure becomes more sluggish in responding 
to a changing situation. One way to overcome this sluggishness is to modify the prior of 
the number of units that are at a particular node. Specifically, suppose there is a sensor 
observation x(n, t + 1) at node n at time t + 1. Before applying the updating procedure of 

(7.5) - (7.7) modify n(k, n, t) as follows. For n > 1 

H{k\n,t) = (l - a(p)) x{ k; n,t) + a[p ) n(k\n - 1,/) (7.10) 

where 0 < a(p ) < 1 is a constant. The constant a(p) may or may not be the same as the 
constant discount factor a(a). If n = 1 (node 1), then 


Jt{k;n,t) = (l - a{p))n{k\n,t) + a(p) s 0 (k ) 


with 


£ 0 (k) = 


1 if* = 0 

0 if it > 0. 


Then, as before 


7 r(k;n,t + l) = C7r(k;n,t)w(k;n,t,x) 
where w(k; n, t, x) is defined in (7.7). 

If there is no sensor observation, one possibility for updating is 


(7.11) 


(7.12) 


n{k\n,t +1) = (l- cto(p))7z(k\n,i) + cc${p)u{k) (7.13) 

where u is a discrete uniform distribution over the possible number of units and 
0 < ao(p) < 1 is a constant. 

b. Modifying the Distribution of the Number of Assets 

The procedure of Section 7.3 updates the conditional mean, m^, and variance, vj , of 
the number of assets at the node given the number of units at the node is k for all 
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possibilities of number of units. It may be advantageous to occasionally reset the 
moments of the number of assets to the TOE (Table of Organization and Equipment) 
values, possibly reduced by perceived attrition. One procedure to reset the moments is 
described in Appendix A. 

Another procedure follows. Associate with each node, not only the posterior mean 
and variance at time 1 for k units, mj^t) and v* (l) but also the TOE mean and variance ky 

and ka 2 a of the number of assets at the node given the number of units that are at the 
node is k. Let 


n 


J(k,n,t) = \ 


0 


if the sensor observation at time 1 belongs to the prior 
distribution of the number of assets for k units at 
node n at time 1 

if the sensor observation at time 1 belongs to the TOE 
distribution of the number of assets for k units at 
node n at time t. 


Let 


ftj{j,k;n,t) = Pjj = j,K{n,t ) = &|all sensor observations during [0,l]j. 

Suppose there is a sensor observation x(n, 1 + 1) =x of the assets at node n at time 
1+1, then a similar procedure to (7.7) can be used to update kj, i.e. 


xj(j,k;n,t + l) = Cnj(j,k;n,t)wj(j,k;n,t,x) 


(7.14) 


with 


wj(j,k;n,t,x ) 


1 

PYn 

1 

s 

1 

l_ 

)) 2 11 

j - OAjJ 

2 

r 2 + v ^ ) [ 

a ( a ) JJ 

1 

1 

1- 

<N 

1 

-\ 

_ I - 

t 2 + kc 2 

2 

[r 2 +ko z a ) 



if .J = 1 

(7.15) 


if / = 0. 
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Further, 


p{k = Ar|all observations in (0,r + l]j = ^7Tj(j,k;n,t + 1) (7.16) 


respectively, 

\ 

P\j = yjall observations in (0,t +1]} = £ Kj(j,k;n,t +1) (7.17) 

k 

gives the posterior probability of there being K = k units at node n at time t +1 
(respectively the posterior probability that the TOE moments summarize the observation 
better than the prior moments). 

Finally, if Pj J = 0|all observations in (0,t + l]J is larger than some constant aj (e.g. 
ay = 0.9) one can reset the posterior moments equal to the TOE values; that is 


m/ c (n,t+ 1 ) = ky 
and 

v 2 k (n,t + \) = k(T 2 a. 

In addition, one can modify the prior distribution of the number of units a la (7.10) - 
(7.13). For example, for n > 1 

Kj{j,k;n,t ) = (l - a(p))7Tj(j,k;n,t) + a{p)rtj(j ,k;n - l,r) 
where 0< cdp) < 1. 


7.5 A Minimal Updating Procedure 

A minimal updating procedure might always use the TOE values for the mean and 
variance of the number of assets and treat the problem of inference of the number of units 
at a node as a classification problem. 
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In particular, let XM(k, n; t) be the estimate of the probability that there are k units at 
node n at time t given all the sensor observations up through time t. Suppose there is a 
sensor observation x(n,t+ 1) of the assets at node n at time r+ 1. Modify the prior 
distribution of the number of units at the node as in (7.10) - (7.11); i.e. 

x M {k,n;t) = (l - a(p))x M (k,n;t) + a(p)?r M (k,n - l;f) (7.18) 

for n > 1 and 


n M {k, 1; t) = (l - a(p))n M (k,l; t ) + a(p)s 0 (k). 

To compute the estimate of the distribution of the number of units at the node 

1 (x(n,t + i)-k/f 


(7.19) 


Ku{k,n-,l +1) = C/r «;r +1)- , 

^r 2 +k<j(af 

where C is a normalization constant. 


exp 


2 r 2 


+ kcr(a) 


(7.20) 


7.6 Examples 

The example of Section 7.2 is simulated. The random times, Tu i = 1,2 of occupancy 
at node i are independent having normal distributions with mean 5 and standard deviation 
1. The standard deviation of the sensor is r. The TOE mean for one unit y— 100 and the 
standard deviation <j a = 10. There is one sensor measurement at node 1 at times 1-7. 
There is one sensor measurement at node 2 at times 5-12. There is one measurement at 
node 3 at times 10-20. The simulation replication is for a total of 20 time units. The true 
number of units is 2. If there are no observations of a node during a reporting cycle, the 
distribution of the number of units at the node is a discrete uniform over 0,1,2,3. 

The following statistics are collected for each replication. If a sensor observation 
occurs at node n, the posterior probability of the true number of units at the node is 

collected. Also collected is the posterior mean number of assets, e.g. 
a = ’^ j mk(n,t)n:{k;n,t). Finally, the average posterior probability of the correct unit is 
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computed over all observations of all nodes over all observation times is computed. Also 
computed is the average absolute deviation between the posterior mean number of assets 
and the true number of assets at the node over all nodes with observations and 
observation times. 

There are 20 replications for each simulation. Table 7.1 displays means and standard 
deviations of the replication statistics over the 20 replications. The results indicate that 
modifying the prior of the number of units as in (7.10) is associated with the greatest 
improvement in the updating procedure for the perception of the number of units at a 
node. If the sensor observation is very good, with small standard error, there is the 
suggestion that resetting the conditional posterior moments of the number of assets to 
their TOE values is worthwhile. Finally, the minimal updating procedure with the 
moments of the distribution of assets always equal to the TOE values appears to be 
adequate. 

7.7 More Than One Observation Per Reporting Cycle 

In the previous subsections of this section there has been 1 observation per reporting 
cycle. In this subsection suggestions for updating procedures with more than one 
observation per reporting cycle are given. 

Since it is known that the ground truth at a node will change, it may not be prudent to 
weight all observations obtained during a cycle equally. Presumably older observations 
should be down-weighted. 

Suppose that observations * 1 , X 2 , ...,xo are taken during a reporting cycle at times t\, 
t 2 > •••> to with t\ < t 2 ^ ... <to- Observation Xj is an observation from a sensor with 
standard deviation % The reporting cycle is assumed to start at time 0 and end at time t. 


60 




a. Summarizing the Sensor Observations 

One possibility for updating is to summarize the O observations and then use the TOE 
classification procedure of Section 7.5 to update the estimate of the probability that there 
are k units at the node at the end of the reporting cycle. 

Table 7.1 


Summary Statistics for Ui 

Ddating Procedures 


Smoothing 

Smoothing 

Posterior 

Average of Mean 

Average of Mean 

in 

Parameter 

Parameter 

Probability 

Posterior 

Deviation of 

■ 

for Number 

for Number 

Used 

Probability of 

Posterior Mean 


of Units 

of Assets 


Correct Number 

Number of Assets 

■ 

dp) 

a(a) 


of Units 

from Correct 

1 

(7.10) 

(7.5) 


(std dev. over 

Number 





simulation repl.) 

(std dev. over 

■ 





simulation repl.) 

20 

0 

1 

n 


41.3 





1 

(5.9) 

20 

0 

1 

nj 

0.72 






(0.06) 

(8.3) 

20 

0 

1 

m 

0.76 

44.3 





(0.07) 

(9.8) 

20 

0.1 

1 

7C 

0.98 

25.7 






(3.8) 

20 

0.1 

1 

XJ 

0.98 

7.37 






(3.4) 

20 

0.1 

1 

m 

0.99 

12.1 






(7.4) 

50 

0.1 

1 

n 

0.82 

28.6 






(6.3) 

50 

0.1 

1 

xj 

0.82 

25.3 






(9.3) 

50 

0.1 

1 

m 

0.81 

22.5 





(0.11) 

(8.7) 

50 

0.1 

0.9 

Tt 

0.79 

28.8 






(7.0) 

50 

0.1 

0.9 

■ xj 

0.79 

25.7 





(0.14) 

(8.8) 

50 

0 

0.9 

n 


49.0 






(6.6) 

50 

0 

0.9 

xj 


54.0 






_(8-7) 
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One weighted average estimate of the number of assets at the node at the end of the 
update cycle is 


x 


a = 


La 2 ( 

/=i Ti a(a) K ' 


o 

i 


l 


(7.21) 


X-i 2 1 \-(‘-‘i) 

/=i rf a(a) y ’ 

This procedure will down-weight older observations. The parameter 0 < a(a) < 1 is a 
smoothing parameter. The smaller a(a) is the less the older observations will influence 
the estimated number of assets, a. 

_2 


Var [A] = 


o 

y— 

'=’ [t? 

72 


o 


(7.22) 


i=i rf a(a) ^ ^ 


Note that if a(a) = 1 and r, = r, then a is the sample average of the observations and 

v 4s]=-d. 

To update the probability of k units at the node at the end of the reporting cycle 


7To(k;n,t) = Cno{k\n,t- A)- 


=T ex P 


1 (a-ky) 1 


2 k(j(af + Var^J 


(7.23) 


V2^-^A:cr(a) 2 + Var^j 

where C is a normalization constant, k 0 is the (possibly modified) probability from the 
previous reporting cycle, and A is the length of the reporting cycle. 

It may be advantageous to have the value of the smoothing parameter aka) depend on 
evidence that the observations, x\, ...,xo are not from the same distribution. If x\, ...,xo 
come from a distribution with the same mean, then 


A+l ~Xi+\ -Xi 


(7.24) 


has mean 0 and variance rj+i + zj 
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Let 1 > ada) > add) > 0 be constants. Set 


a(a ) = 


a m (a) if 


a c {a) if 




7+1 


yR+i+rl 


■<c(c) 


lA'+i 


(7.25) 




_2 - 2 
*i+l ^z 


>c(c) 


where c(c) is a constant; for example c(c) = 3. 


7.8 Examples 

The simulation model of Section 7.6 is run. Observations occur at unit times. All 
observations have sensor standard deviation r. The reporting cycle time is either 2 or 5 
time units. The parameter indicator of change, c(c) and the asset smoothing parameter if 
change is indicated, add) are varied. The TOE classification procedure (7.23) with n 0 
modified as in (7.10) - (7.11) is used to estimate the number of units at the node. The 
total number of replications is 100. The summary statistics of Section 7.6 are gathered. 
The results appear in Table 7.2. 

The results suggest that summarizing the data with an average of the observations is 
not as good as a weighted average with older observations having less weight. The results 
also suggest that using information from the previous node at the end of the previous 
reporting cycle is advantageous. 

The physical movement and observation model of Section 7.6 is simulated again. 
However, there are two different sensors; one with standard deviation zi = 20 and the 
other with standard deviation th= 50. Observations are taken at integer times and the 
sensor used is chosen at random. Each sensor has probability 1/2 of being chosen. If two 
observations are taken at a time both observations use the same sensor. 
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Table 7.2 

Summary Statistics for Updating Procedures 
With More Than One Observation Per Reporting Cycle 



Reporting 

Cycle 

Time 

Indicator 

of 

Change 

C(C) 

(7.25) 

Asset 
Smoothing 
Parameter 
No Change 
ania) 
(7.25) 

Asset 

Smoothing 

Parameter 

Change 

<*d°) 

(7.25) 

Smooth 

Parameter 

for 

Number of 
Units 
<dp) 
(7.10) 

Posterior 

Probability 

Used 

Average of 
Mean 
Posterior 
Probability 
of Correct 
Number of 
Units 
(std dev. 
over 

simulation 

repl.) 

Average of 
Mean 

Deviation of 
Posterior 
Mean 
Number of 
Assets from 
Correct 
Number 
(std dev. 
over 

simulation 

repl.) 

20 

2 

3 

1 

0.1 

0.1 

no 

0.87 

(0.11) 

32.7 

(13.4) 

20 

2 

3 

1 

0.1 

0 

no 

0.57 

(0.11) 

67.6 

(11.7) 

20 

2 

3 

1 

0.1 

0.2 

no 

0.85 

(0.12) 

33.5 

(13.1) 

20 

2 

3 

1 

0.05 

0.2 

no 

MB 


20 

5 

3 

1 

0.05 

0.1 

no 

0.88 

(0.12) 

46.7 

(24.2) 

20 

5 

3 

1 

0.10 

0.1 

no 

0.86 

(0.13) 

48.6 

(20.2) 

20 

•5 

2 

1 

0.05 

0.1 

no 

0.85 

(0.12) 

52.1 

(21.1) 

20 

5 

3 

1 

1 

0.1 

no 

0.70 

(0.18) 

71.0 

(19.9) 

20 

5 

3 • 

0.05 

0.05 

0.1 

no 

0.87 

(0.12) 

48.2 

(22.6) 

50 

5 

3 

1 

0.05 

0.1 

no 

111 

63.2 

(24.3) 

50 

5 

3 

1 

0.05 

0 

no 

0.46 

(0.15) 

80.8 

(27.8) 

50 

5 

3 

1 

0.05 

0.2 

no 

0.63 

(0.12) 

64.0 

(21.5) 

50 

5 

3 

1 

0.1 

0.1 

no 

0.62 

(0.13) 

66.8 

(22.0) 

50 

5 

2 

1 

1 

0.1 

no 

0.62 

(0.14) 

74.5 

(18.3) 

50 

5 

3 

0.05 

0.05 

0.1 

no 

0.60 

(0.15) 

65.8 

(24.9) 
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Table 7.3 reports the simulation results. Each simulation has 100 replications. Once 
again, there is the indication that weighted averages of observations with older 
observations weighted less are better than simple averages. There is also an indication 
that using prior information from neighboring nodes is better than only using information 
from the node. The values of the various smoothing parameters do not seem to be that 
important. 


Table 7.3 

Summary Statistics for Updating Procedures 
_Randomly Chosen Sensors_ 



Asset Estimator Parameters 

Unit 

Estimator 

Parameter 



Reporting 

Cycle 

Time 

Indicator 
of Change 

Asset 
Smoothing 
Parameter 
No Change 
aN(a) 

Asset 

Smoothing 

Parameter 

Change 

Smooth 

Parameter 

for 

Number of 
Units 
a(p) 

Mean 

Average of 
Posterior 
Probability 
of Correct 
Number of 
Units 
(std dev. 
over 

simulation 

repl.) 

Mean Average 
of Deviation of 
Posterior Mean 
Number of 
Assets from 
Correct Number 
of Assets 
(std dev. over 
simulation 
repl.) 

5. 

3 

1 

0.1 

0.1 


58.1 

(22.7) 

5 

3 

1 

0.1 

0 

0.55 

(0.15) 

88.5 

(24.2) 

5 

3 

1 

1 

0 

0.47 

(0.12) 

97.1 

(19.4) 

5 

3 

1 

1 

0.1 

0.70 

(0.15) 

71.9 

(21.0) 

5 

3 

0.1 

0.1 

0 

0.57 

(0.19) 

77.3 

(31.3) 

5 

3 

0.1 

0.1 

0.1 

0.76 

(0.13) 

55.2 

(21.7) 
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7.9 Generalization to More Than 1 Asset Type 

Suppose that a unit has / asset types. Let Xjj denote the z 1 * 1 sensor observation of asset 
type j. Let rji be the standard deviation of the z* sensor observation of asset type j. Let tp 
be the time of the sensor observation and Oj be the number of sensor observations of the 
7 th asset type during the reporting cycle. A weighted average estimate of the number of 
assets of type j is 


Oi 


K Jl 


a / = 


and 


Va i[Aj] = 


1=1 Tjia(a) 


Oj 1 

y 1 


i=i r)ia(a) 


Oj 

v 

r 2 ji 

/ 

1=1 

a)-^f 

(z 

• T 


(7.26) 


(7.27) 


2 / \-U~tji) 

i=i Tjia{a) v 1 ’ 

To update the probability of k units at the node at the end of the reporting cycle 
7u 0 (k;n,t) = Cn 0 {k-n,t- A)J| , = = exp{-± - r -— \ ( 7 - 2 8) 


l 


rexp< 

1 (dj-kyjf 


ij kcr a{j) 2 + Var 

A 

2 ka a (jf + Var 



where yj is the TOE mean number of assets of type j for 1 unit, cr a (j) is the TOE standard 
deviation of the number of assets of type j for 1 unit, C is a normalizing constant, 7co is 
the (possibly modified) probability from the previous reporting cycle, and A is the length 
of the reporting cycle. 
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7.10 Updating Using Same Node Information 

Consider just one asset type at a node and let /u(t) be the ground truth number of assets 
of that type on the node at time t. Assume observations of the asset number are taken at 
times t\ < t 2 < .. .<t n during an observation period. Let X( 4 ) denote the observation 
obtained at time 4 . Assume 

X(*k) = /<4) + ®(4) 

where fE( 4 ) has a normal distribution with mean 0 and variance z( 4 ) 2 ; the observation 
X( 4 ) could also be obtained from the binomial or multinomial sensor models of Sections 
4 - 6 ; the appropriate variances should be used in the procedures below. 

Estimates of node occupancy are obtained at the end of observation periods and use 
sensor observations collected during the period. Since units move, ground truth /u(t) will 
change as a function of t. The problem of detecting when changes occur in ground truth 
from observational data has received attention for many years; cf. Basseville (1988) and 
Lai (1995). Problems of this type are called changepoint problems. If a change in ground 
truth occurs during an observation period, then combining observations as proposed in 
Appendix A may not be the best procedure. 

There will be many nodes and arcs in a typical theater-level model. Sensor 
observations may occur at many nodes and more than once at a node during an 
observation period. As a result, the statistical calculation to estimate the number of assets 
at a node will be done many times for each observation period. Thus, one must be careful 
with regard to the computational burden of the statistical calculation. Another procedure 
for combining observations is proposed below. 
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A proposed procedure for combining observations which may be from distributions 
having different means 

We will take a generic observation period to be one time unit in length; e.g. [0,1]. Let 
c be a tuning parameter chosen by the analyst; good values for c are 2 and 3. 

Suppose observations are taken at times 0 < t\ < t 2 < ...<t n < 1. 

1. Compute 


n jjfrMjM . 

, 1z( 4-,)-X(/ b _ 2 )1 


D J = 




)+T 


(0-0 


2. If max (£),•) < c, then compute 

j = n ,...,2 


(7.29) 


y 4?*) 

-.iiwi 

* *v*) 


and v = 


y_L_ 

rr(4) 2 


(7.30) 


Use /i as the estimate of the number of assets at the node and v as its 


variance. 

3. If max (Dj) >c, then find the largest j for which D, > c. 

J- max{/: Dj > c}. 

Compute 


y 44) 
. *>y^(4) 2 

v= —y- 

Sr(4) 2 


and 


v = 


1 


E 

Ax/ 



(7.31) 


(7.32) 
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7.11 Procedures for Updating Number and Type of Unit at a Node 

For simplicity of discussion and notation we will assume there is one type of unit with 
only one type of asset. Let K(n, t) be the number of units that occupy node n at time t and 
let A(n, t) be the number of assets at node n at time t; n could also be an arc. Assume the 
prior distribution for the number of assets given the number of units with no sensor 
observations is 


P{A(n,0)eda\K(nS)) = k) = 


J2^4ko a H 


J .{a- kyf 
2 kcrl 


(7.33) 


where y is the TOE (Table of Organization and Equipment) value for the amount of asset 
for one unit. 

Let ’ 


7t(k;n,t ) = P^K{n,t) = &| all sensor observations during [0,r]} (7.34) 

and let 

dF k (a;n,t ) = p{A(n,t) eda\ all sensor observations during [0 ,t],K(n,t) = 

1 f 1 ' 2 ) (7-35) 

"&v 1 (»,,)T2 (a '"‘ (M)) lvl(n ' t) \ 

Results from a small-scale simulation reported in Section 7.6 suggest that modifying 
the prior distribution used in the Bayesian updating of the distribution of the number of 
units at the node can improve the performance of the updating procedure. The 
modification uses information concerning the posterior distribution of the number and 
types of units at the nodes which are direct neighbors of n; the specification of direct 
neighbors can depend on the length of time between sensor observations. 

Suppose that node n has only 2 neighbors «_ and n+; these neighbors could be 
adjacent arcs. Modify the prior distribution of the number and type of units at node n at 
the beginning of the Z* observation cycle as follows. 
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GC CC CL CL 

n{k,n,t ) = — n{k,n-,t) + (l - a)n{k,n,t ) + — n{k\n+,t) + — f 0 (£) +—/?(&) ( 7 . 36 ) 

4 4 4 4 

where a is a tuning parameter; might try a = 0.1; so(k) = 0 if k * 0 and £o(0) = 1; and j6 

assigns equal probability to each possible number and type of unit combination. 

A 

Let A be the estimate of the number of assets at node at the end of the observation 

A A 

cycle; let V be the variance of A . The posterior distribution of the number and type of 
units at node n at the end of the observation cycle is 


+ l) = Cx(k,n,t + l) 


1 


2 KtJV + kcr(a ) 


ex P -r-7 


(A-ky) 


2 1 


2 F + kcy(a ) 


(7.37) 


notice that the TOE values are used rather than a current estimate of the mean and 
variance of the number of assets at the node given the number of units. 

Numerical computation considerations suggest replacing (7.37) by the following 
computation. 

1. For each k, compute 

f{k) = \oz7r(k,n,t + \)-\\og{y + M q ) 2 )- "| ~ , yL • 

z V + ka\a) 


2. Find 


3. Let 


f m = max f(k). 

k 


K = {k-\f(k)-f m \<M} 


the set of all number and types of units whose / value is within M of the 
maximum; the value of M can be chosen by the analyst; M-4 might be 
reasonable. 

4. For keK, let 

n(k,n,t + \) = 


exp{/(£)} 

X ex p{/1/)} ‘ 

i — v 
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For k&K 


n{k,n,t + l) = 0. 


7.12 Updating Using Same Node Information with Binomial-like Sensors 

Consider just one asset type at a node. Let /u(t) be the ground truth number of assets of 
that type on the node at time t. Assume observations of the asset number are taken at 
times t\<t2<...<tK during an observation period. Let X(tk) denote the observation 
obtained at time tk- Assume the conditional distribution of X(tk) given ju(tk) is binomial 
with ju(tk) trials and probability of success p(tk). In Section 4, an estimate of p(tk) is 
proposed in (4.19) 

,_X{t k 

k) ~ pW)' 

It is suggested that the distribution of ju(t k ) be approximated by a normal distribution 
with variance r 2 (tk) = X(tk)( 1 - p(tk))/p(tk)'> if X(t k ) = 0, then set r(tk) 2 = vo > 0 where 
vo is chosen by the analyst. A modification to the procedure proposed in Section 7.10 is 
as follows. Let c be a tuning parameter chosen by the analyst. 

1. Compute 

X(t k ) X{t h - 1) 

D P^k) P{h~\) 

^{tk) + r 2 (t k .i) 

X(t k . i) X(t k - 2 ) 

D k pjtk-l) 
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2. If max (D y ) < c, then compute 

j~K,...y 2 


Y x{tk ) 1 

k p(tk) r(t k ) 2 

Zih 

k T(t k ) 


and 


v = • 


1 ' 


7r(4) 

Use ji as the estimates of the number of assets at the node and v as its 


variance. 


3. If max (Dj) > c , then find the largest j for which Dj > c. 

j=K,..., 2 


Compute 


J=max.{j: Dj >c}. 

y *(4) 1 

k>jp{ { k) r(t k f 


M = - 


V 1 

fcjr{tk) 2 


and 


v = 


^ ) 2 


8. Ground Course-Of-Action Realization and Perception 


8.1 The Problem 

A theater-level action, in simplest form, involves the movement of units of a force, 
here called Blue (B), such as infantry or armor brigades or divisions, from their initial 
locations (network nodes) to a final destination (possibly a single network node). The 
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ground forces may be supported by air reconnaissance and strikes. The course of action 
(COA) plan could be that the units arrive at the destination, D, at a particular time and in 
coordination, although quite possibly along different routes or paths. In reality, times of 
transit from place to place (node to node) may vary, e.g. for reasons of weather or natural 
terrain obstacles or equipment breakdown, etc., so that times of transit from node to node 
may vary randomly. The arrival times at, and directions of approach to, the destination 
may vary as well, but by design of B-force commanders. Of course the rate of advance 
and force allocations depend upon B perceptions concerning the opposition’s actions and 
capabilities. 

From the point of view of the opposition, here called Red (R), the options for a B 
action are known broadly. As the campaign develops the actual COA being carried out by 
B will become increasingly evident to R. Of course R has his own intended COA, which 
may develop spatially in conformity with what his C3I system tells him about B. 

Given the above general background we wish to provide an analyst with the capability 
of obtaining and maintaining a quantitative perception of one opponent’s COA (here B) 
by the other (here R), based on the latter’s sensor information and prior probabilities of 
the various COAs. The methodology applies to either side. 

8.2. Routes and Corridors of Advance as Part of a Course of Action 

The definition of a COA involves a mobility corridor, which may comprise several 
specific routes, e.g. highways or trails and neighboring territory over a broad geographical 
area; it may even involve a sea-land amphibious landing action. Only one route may be 
used by all units in the particular COA, or several routes may be used by different 
numbers of units. It is assumed that elements of major forces such as brigades or 
divisions (units) move on such defined routes, i.e. along one or more axes of advance at a 
speed influenced by the type of unit and the nature of the route. It is a convenient 
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modeling assumption that units are detected on such routes; the probability of detection at 
any time along a route (an arc, or node, in the current model) depends on the unit sizes 
and types (their asset portfolio) occupying the route at the time and the types of 
surveillance assets (sensors) viewing the route at the time. 

Consider the problem of detecting and classifying an advancing force. It seems 
reasonable to think in terms of sensor updating at regular time intervals A apart, where 
A = 6 to 12 hours when thinking about ground unit advance. Suppose no units have been 
detected up to time kA on the routes within a given corridor. The observing force must 
decide where to place his reconnaissance effort within the particular corridor for the 
following time period (kA,(k + 1) A], and he must update his probability of corridor 
usage depending upon reports received from his reconnaissance system within that time 
period. One way of proceeding is to specify that particular units will use particular routes 
within the corridor, advancing at specific speeds. If they start the journey at specified 
times then it can be forecast where they will be at specified times in the future. This 
implies that times of exposure on particular parts of a particular route during time interval 
(kA, (k + 1) A] are specified. If reconnaissance assets are assigned to those locations 
during those times the only reason that the units will not be detected is through failure of 
the assets to perform adequately, given the presence of the units sought. 

Uncertainties 

Inherent in the above situation are various realistic variabilities and uncertainties. 
First, the exposure times of particular parts of a particular route (arc segments and node 
residence times in the arc-node model) will vary because start times and transit times 
along route segments (arcs, nodes) will tend to vary, and are actually modeled as varying 
randomly. Note that the probability law of transit and residence time variation that 
governs movement is supposed known to the agent performing surveillance; this is 
optimistic and can be changed. It may be that the maneuvering force will deliberately 
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increase the variability of its movement so as to throw off the enemy reconnaissance 
effort. 

Second, uncertainty exists concerning reconnaissance or surveillance asset 
effectiveness: successful detection performance is not certain, even if assets are employed 
during a maneuvering unit’s exposure time in a particular geographical region. The 
probability that the reconnaissance asset flies close enough to the force to detect etc., is 
not unity; the probability that overhead surveillance assets detect may not be unity 
because of cloud cover or weather conditions. Realistically, too, delays will exist in 
processing detection reports and in corroborating them. 

A third uncertainty factor exists concerning the size and asset composition (e.g. tanks, 
APCs, infantry) of the maneuvering units. 

8.3 Analyst Tasks and View 

In what follows we describe the way in which an analyst must proceed to use 
J-STOCHWARS; we attempt to describe this stepwise, although there probably are gaps 
that must be filled. In particular we focus on B’s ground force movement, and on R’s 
perception of it by reconnaissance assets that are mainly airborne or in space, until Blue’s 
objective is reached, where ground forces are assumed to be present. 

Analyst Steps: Modeling Red’s Perception of Blue’s CO A. 

(1) Establish Spatial Network: physical nodes and connecting arcs 

(2) Specify current objective for (e.g.) Blue and a schedule by which forces 
(specify types) move from Source nodes to (current) Destination. 

Also for opposition, Red. 
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D (Red) 



(3) Initial/Blue COA Determination: analyst specifies geography of the advance, 
options, and the numbers and timing of units advancing, first without 
opposition (certain forms of opposition can probably be optionally 
automatically represented, e.g. mine fields or even airstrikes). Broad mobility 
corridors will be set up for B, within which more specific axes of advance are 
specified. Red (R) must, in turn, assess B’s avenues of approach. The model 
first takes these to be the same as the axes of advance. The B units will 
advance along arc-node paths within these. 

(a) Specify nodes & arcs in each axis of advance (see above). (Tentative; analyst 
intervention). 

(b) Determine spatial realizations of each a.a. (avenue of approach): define 
possible routes for forces within bands defining mobility corridors. Develop 
alternatives within bands. 

(b-1) Routes and force types are not independent: some routes are appropriate 
only for certain force types. 

(b-2) Optional routes, and forces thereon, could be specified by analyst. 

Consider location of opposition in picking routes and timing advance. 
(b-3) Optional routes can possibly be specified “automatically”, e.g., by 
randomization with probabilities specified by the analyst over particular 
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route options. The analyst, acting for Red, must first develop a portfolio 
of routes with possible unit types on them. These are Blue COA 
options. 

(c) (Blue) picks one COA option that she will actually use; i.e., a bundle of 
paths/routes from each movement band consistent with a COA. These can be 
called sub-COAs. In general there can be several sub-COAs, i.e. alternative 
routes and forces thereon within a given COA. 

(c-1) Respects route-asset compatibility constraints. 

(c-2) Done by analyst intervention or by randomization. 

(c-3) Blue remains on route until he encounters, or anticipates, opposition. At 
some point may change COA. The present discussion does not yet 
include encounters with opposition. 

(d) All of above is done on the basis of initial perception of opposition’s (Red’s) 
locations and force composition and strength. It is subject to change. Success 
depends on the quality/accuracy of that perception. 

(4) The Blue COA option is now simulated: transit times on arcs are simulated 
from arc-time distributions. 

(a) Detection times by R sensor assets on arcs are simulated, potentially providing 
information on location and type of B units and assets thereof to R; detections 
occur at an exponential rate during units’ exposure time to the sensors; the 
exponential rate will be a function of sensor allocation, detection range, etc.; 
time-dependence is allowed by letting detection rate on an arc depend on the 
time. 

(5) Next develop quantitative (Red) perception of the opponent’s (B’s) COA. This 
develops from sensor data and prior information. 

(a) (Red) allocates sensor effort to arcs/nodes in (Blue) COA’s. Options : 

(a-1) Analyst specifies (using experience), with the aid of the Air Model. 

Translate into detection rates of various units/assets on various arcs. 

(a-2) Automatically determined using an optimization model (an option for 
the future). 
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8.4 A Proposed COA Perception Update Procedure Based on Unit 
Detection. 

In this section we detail an approach to COA quantitative perception updating. The 
approach involves simulation in order to avoid practically infeasible tailor-made 
analytical calculations. We first describe procedures that update perception based on 
whether or not units are detected. Analytical calculations are computationally difficult 
because COAs are not specified by single paths, but rather by mobility corridors which 
can contain many possible paths. For a realistic network, the number of different ways 
units may traverse a network for one course of action will be large. Furthermore, paths 
belonging to different COAs may overlap during parts of the scenarios. Hence an analytic 
approach to updating perception based on enumeration of all such paths does not appear 
practical. 

(1) Initial Steps. 

(a) The analyst determines a COA perception update cycle time, A. 

(b) At times kA(k — 1,2,...) measured from an initial time (e.g. action beginning) 
the probabilities of the opponent’s (here B’s) COAs are calculated; those at 
time (k + 1)A are obtained from those at time kA, augmented by information 
obtained in the time interval (kA, (k + 1)A]. In practice A might be 6-12 hours 
in duration. 

(c) Let Q be the event that the z* (i = 1,2,..., 7) COA is being pursued (by B); let 
Q(k) be the event that Q is perceived (here by R) to be in progress at time k. 

Then let n,(^) - , the posterior distribution for various COAs at 

time kA. This quantifies R’s current (time k ) perception of B’s forces and 
COA. 
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D (Red) 

© 



C 2 1 0 5 


Figure 8.2 

(d) Consider a given COA, so let Q begin at time 0. Suppose for illustration we 
have the two alternatives shown in Figure 8.2 and that C\ is actually being 
pursued. Ci has 4 units traveling along the leftmost route in Avenue 1, 1 unit 
traveling along the rightmost route in Avenue 1 and 1 unit traveling along 
Avenue 2. At time Jne[0, A) a detection occurs on the left-most route in 
Avenue 1 (denoted An); at time e?i2e[0, A) a detection is made on An; no 
detections are made on^2l> the only path in Aj. 

The above corresponds to actual data obtained in a real operation. The 
objective is to obtain from it a probability that C\ is actually under way, plus 
probabilities of alternatives (here only C 2 ). 

(2) Evaluating the Likelihood of a COA. 

We want to evaluate the probability of observing what has been observed in terms of 
the possible COAs, i.e. the likelihood that one or the other COA is in force. This can then 
be combined with TLj(k) to update to n z (A + 1). We want to do so economically, so we 
choose to approximate and simplify by removing the specificity of the observations 
mentioned in (l),(d): we report seeing 2 detections in Avenue 1 and 0 in Avenue 2. 
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We now calculate the likelihood of Cj and Cj on the basis of the above data 


presentation. Here are two procedures. An exposure time during (kA, (k + 1)A] is the 
amount of time in (kA, (k + 1)A] that units on an arc are subject to sensor observation; it 
is a function of transit times, the type of sensor, and environmental variables. Let k\\(k) 
be the exponential rate of detection for one unit during its exposure time on A\\ in 
(kA, (k+ 1)A]. The following is a simulation approach to evaluating the likelihood of a 


COA. 


(a) 


Simulation 

(a-1) Suppose Cj is in effect. Simulate the exposure times, on each arc, of 
4 units on A u and 1 unit on A \2 during (0, A], and also simulate for 
each exposure time the number of times detections occur. Note that the 
detection rate on An is X\\(k)4 and on A \2 is Xu(k)\ for Cj; it is 
A\\(k)\ forC2. 

Use the constraint that if a detection occurs on an arc in an avenue of 
approach no more can occur thereafter on that arc or on routes/paths 
that continue it (unless the arc ends at a node out of which several arcs 
occur; in that case the units originally detected could be lost, and might 
need to be re-detected). There can thus be 

«00 replications with no detections on either An orAu; 
n 10 replications with 1 on A 11 , 0 on A 12 ; 
no 1 replications with 0 on A 11 ,1 on A \ 2 ; 
n 11 replications with 1 on A \ 1 ,1 on A 12 ; 

(a-2) To estimate the likelihood of Cj given 2 detections on A\ (one on A\\ 


and one on A 12 ) quote 


»n(l) 
. «' 0 ) 


V 


; here «n(l) refers to one detection on 


each path for A\ and « ; (1) is the number of simulation replications. 
Presumably one also simulates events on A 2 , using A.2'1 for the 

«o(2) 


detection rate for Cj; the probability of the observed data is 


n,(2) 


where «o(2) is the number of n^2) replications that result in 0 detections 
on A 2 . 
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Hence the (combined, estimated) likelihood for C\ is cdC \; data) = 

\i(i) Y ”0(2) ^ 

(a-3) Next re-simulate the detections but with the detection rates appropriate 
to COA C 2 

C 2 : /lirl on(and A 12 ); 

A 2 -5 on A 2 

and tabulate the same simulated ratios. Apply Bayes to update. 

(a-4) The above procedure will work for arbitrarily complicated setups, but 
may be computer intensive. An alternative is 
(b) Hybrid Simulation and Calculation of Likelihood. 

It is advantageous to avoid the simulation of detections in the process of 
updating perceptions of alternative COA probabilities. A possible way of 
doing so is as follows. 

(b-1) Consider the r* exposure time realization for arcs and nodes occupied 
during (kA, (k+ 1)A]. Let 

M k ’ r )= T i A j( k X k > r )n u A k )/ A (s.i) 

jzj(k) 

be the mean detection rate over the entire mobility corridor for period k : 
times (M, (k+ 1)A]. In (8.1) we sum up the basic detection rates on all 
arcs and nodes that are occupied (have non-zero exposure times 
e(k;r)ji) during period k)\ these are weighted by the appropriate 
exposure times for COA Q, and by the arc/node loads (unit sizes), 
Uji(k) that appear on arc/node j for COA Q. 

We now compute that the probability of 0 detections as 
P{D(k) = 0\Q,r} = e~ Ji{k;r)A 

This is the likelihood element for C ; - given the exposure times of 
realization r\ if exposure times are in effect sampled independently we 
use 

K r=1 

where R is the number of replications. 
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(b-2) Suppose at least one detection takes place, an event of probability 
1 — e _A '^*’ r ) A =Xi(k\r) A if /l,-A is small, as may be true in practice. 

Simple Approximate Approach When Detections Occur in 

(M, (k+ 1)A]. 

Assume that if 1 or more detections occur in the interval 

P{m - ■ 

If the above is averaged over the R replications of exposure time we 
arrive at an expression for 

P{D(k) = d(qQ}sp d{k) (i,k) 

an average of Poisson probabilities. This is the semi-parametric version 
of the Simulation approach above, in (2), (a). 

A More Correct (but Difficult) Approach 

The above procedure is not literally correct (although it may be an 
adequate approximation) since once a unit detection occurs on a 
specific route there are assumed to be no more detections on that route 
unless the unit reaches a physical node where she can be hidden or 
emerge “branched” or “split” into two or more unit segments (e.g. a 
division may divide into several brigades) which take different routes. It 
becomes necessary for the analyst to decide how the split occurs; since 
this can be done in several ways a multiplicity problem threatens; one 
way to handle this is by sampling. A path on which a detection has 
occurred during period k is essentially pruned for the purpose of 
detection after the detection on it in period k ; this holds into period 
k+ 1, or until a node is reached with a split. Note: another way of 
pruning is to ignore all detection events (counts) greater than 1. 

After detection occurs classification assets are brought into play. 
The updating that follows will be described subsequently. 

8.5 Updating the Perception of COA Using Asset-Counting Sensors. 

Once units are detected, asset counting sensors can be employed to obtain more 
information. Suppose units have been detected at node N„. Let Si(n,j; k), ..., St, n ( nj ; k) 
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be the sensor observations of assets of type j at node N n during the observation period 
(kA, (k + 1)A]. We will assume A is small enough so that the number of assets at the node 
is constant over the period. The least-squares estimate of the number of assets of type j at 
node N n based on the latest sensor information only (gathered during period kA, (k + 1)A 
is 




where z 2 nj is the variance of the error of the sensor observation S e (n,j; k). 


We can interpret the results of the sensor observations as resulting in a normal 
distribution for the number of assets of type j at the node during the time period; the 
distribution has mean 


Un / On 1 

m M = / E-277: 

_t=l J/ t- 1 T nj\") 


and variance 


v n j{k)~ ^ - (8.4) 

§ * 2 nj(/) 

Let J(cc, k ) be the collection of all nodes that might be occupied during the period 
(kA, (k+ 1)A] for a particular avenue of approach a. The result of all sensor observations 
during (kA, (k + 1)A] is a distribution of the total number of assets of type j at all the 
nodes in J(oc, k) which is normal with mean 

m.j(a;k)= J^m^k) (8.5) 

N n eJ(a\k) 


and variance 


o 2 j(a;k)= ^o 2 nj (k); 

N„eJ{cr,k) 


( 8 . 6 ) 
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how well this distribution reflects the true number of assets that are on the avenue of 
approach a will be affected by the number of sensor observations at each relevant node. 

A particular COA c has a distribution of the total number of assets of type j at nodes 
N n eJ(or, k) in the avenue of approach a during the time period (kA, (k + 1)A], Aj ( a;k ), 

P{Aj{a\k) eda\COA = = %{a\ fij{a,k',c),G 2 j{a,k-,c))da (8.7) 

where <£( a; pi, o 2 ) denotes the normal density function with mean pi and variance o 2 ; the 
mean pij(a, k; c ) is obtained from the TOEs (Table of Organization and Equipment) for 
the types and numbers of units using the avenue of approach a for the COA c. 

The procedure to update Red’s perception of Blue’s COA is as follows. Let n(c; k) be 
Red’s posterior probability at time kA that Blue is following COA c. To obtain a posterior 
probability which incorporates the sensor information obtained in (kA, (k + 1)A] 

n(c; k +1) = DU(c; Yl J %(a, pij {a, k; c), a) (a; k; cj)%(a, m j (a; k), u) (a; k))da 

« j 

( 8 . 8 ) 

= DU(c, £) J} 11 ( a > k )> t* J ( a ’ k ’ c ')’ v 2 j(a>k) + <j) (< a\k,c )) 

a j 

where the product a is over the avenues of approach, the product j is over the asset types 
and D is a normalization constant. 

Note: the above can be adjusted for the time that the sensors have to “measure” the unit 
assets during (kA, (k+ 1)A]; it can also be made to incorporate a likelihood component 
informative of assets from the detection. Note also that since the procedure is based on 
the total numbers of assets each COA has on different avenues of approach, the procedure 
will produce more discrimination between COAs the more differentiated the COAs are. 
The COAs will be better differentiated the greater the difference of the numbers of assets 
on different avenues of approach for different COAs. The COAs will also be better 
differentiated in situations in which the avenues of approach do not have nodes and arcs 
in common. 


84 



8.6 Numerical Considerations 

The proposed updating of the course of action of (8.8) is susceptible to numerical 
problems. We propose to first compute the logarithm of the normal density functions. 
First compute 

t x (c, k + 1) = £- ]- (nij (a; k ) - (a; k, c)) 2 /[oj («; k )+Oj («; k, c)] (8.9) 

a j 1 

£ 2 {c,k + l) = Y J Tj-\^fij( a > k ) + oj(a;k,c)) (8.10) 

a j 2 


t{c, k + 1) = i\ (c, k + 1) + i 2 (c, k + 1) + In fl(c, k). (8-11) 

Order £(c, k + 1) from largest to smallest; e.g. 

(& +1) = max l{c, A: +1); (8.12) 

C 

call these £(i)(k+ 1), i( 2 )(k + 1),.... Let c(i), C( 2 ) be the corresponding Course of Action 
(COA), e.g. 

C(j) = jc:^(c,A: + l) = max.£(c,A: + l)J; (8.13) 

if there is more than one COA in the set, order them in some fashion. 

Consider those COA’s such that 


+ 1 )-.£(;)(£+ !) 


<4; 


(8.14) 


call these c(i), C(2), ..., C( m y Let R be the number of COA’s not among these m. 

Let 0 < a < 1 be a parameter to be determined by the analyst; e.g. a = 0.9. 
For c€{c (1) ,c (2) ,...,c (m) },put 


fl(c; k +1) = a 


For cg{c (1) ,c (2 ),...,c (m) },put 


exp{^(c, & +1) - (k + 1) j 

X exp{^(y, & +1) - £{\){k +1)} 


(8.15) 


85 


(8.16) 


n(c; £ +1) = (1 - a)—. 

Another possible procedure to update the probabilities of the different COA’s is the 
following. Modify the prior 




(8.17) 


where |Cj is the number of different COA’s and a is a parameter chosen by the analyst, 
e.g., a =0.1. Compute 


£(c, k +1) = £ , (c, k + 1) +£ 2 {c,k + 1) + In fi(c, k). 

Order £(c, k + 1) from largest to smallest and define £(j)(k +1), cq), as before. Let 

, exp{f(c,i-H)-< (1) (* + l)} 

+ 2 « P «/ > t + i)- <w (t + i)} 


for ce{c (1 ),C( 2 ),...,c (m )};for c?{c(,),...,c w ) let n(c, k+ 1) = 0. 


(8.18) 


(8.19) 


8.7 Degrading the Estimate of the Number of Assets on a Node for Which 
There Are No Recent Sensor Observations 

Suppose the last sensor observation of node n occurs during observation period k. Let 
Aj{k\n ) be the estimate of the number of assets of type j on node n at the end of 

observation period k; let o 2 j(k;n) be the variance of the estimate. Suppose there are no 
sensor observations of the assets at node n for the next m observation periods. Since the 
node has not been observed for m observation periods, there may be more uncertainty 
concerning the estimate of the number of assets of type j at node n. Thus we will inflate 
the variance, by setting 

vj(k+m-,n) = ^A ( 8 . 20 ) 

a 
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where 0 < a< 1 is a parameter determined by the analyst. If a= 1, there is no inflation. 
The standard deviation is inflated by (1 /4a ) m . 


Inflation Factor for Standard Deviation of the Estimate: 1/ 4a 


warn 

1 

2 

3 

... 

OO 

... 

16 

0.8 

1.25 

1.6 

1.95 


6.0 


35.5 

0.9 

1.1 

1.2 

1.4 


2.3 


5.4 

0.95 

1.05 

1.1 

1.2 


1.5 




Suppose the observation intervals are 3 hours long. If 4a = 0.9 and there have been no 
observations for 24 hours (8 observation periods), the standard deviation of the estimate 
is multiplied by 2.3. It may be appropriate for nodes and arcs to have individual a’s; one 
a for a node/arc which units travel through; another a for a node on which assets are 
known to stay for a period of time. 

To improve the numerical stability of the COA updating calculation, it may be 
advantageous to first do a test of hypothesis that the estimated number of assets of type j 
at node n is different from zero. 

Let J(cc, k) be the collection of all nodes that might be occupied during the 
observation period [kA, (&+l)A] for a particular avenue of approach. For each node 
neJ(cr, k) and type of asset j, compute 


Oj(k;n) 

where vj(k; ri) has been multiplied by (1/ 4a ) m if node n has not been observed by a 
sensor during the last m observation periods. If 


Bj(k;n)< fi 
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then do not include node n and asset j in the calculation for COA update where (3 is 

chosen by the analyst; (3=2 may be a good choice. If all asset estimates are 0 for J(cc, k), 
let oj be the smallest variance. 

8.8 A Procedure to Update COA, Using Asset Counting Sensors, Which 
Recognizes the Number and Type of Unit on a Node/Arc Changes Over 
the Updating Period 

The estimates of the number and type of units on each node and the number of assets 
of each type on a node is updated at the end of each sensor update cycle; nominally 
2 hours, cf. Sections 7.6 and 7.11. 

The COA update cycle length is often longer; nominally 6-12 hours. 

It is proposed that the total number of assets of type j in the avenue of approach a, 
Aj{a;k ) be estimated as the sum of all the estimates of the number of assets of type j on 

all the nodes and arcs in avenue a obtained in the latest sensor update cycle. The 
estimation procedure for the number of assets of type j on a node/arc would incorporate 
procedures in Section 7.10 that address a changing ground truth during the sensor update 
cycle and that of Section 8.7 which degrades the estimate on a node/arc for which there 
are no recent sensor observations. 

Let m„j be the estimate of the number of assets of type j, on node/arc n at the end of 
the latest sensor update cycle and be its variance. The estimate of the total number of 

assets of type j in a particular avenue of approach a is obtained from (8.5). Finally the 
probability that each COA is being executed is obtained from (8.9) - (8.12) and (8.17) - 
(8.19). 
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TIMES TO TRANSIT ARC 
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Figure 1 

TIMES TO TRANSIT ARC 
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MATMOV[;1] 


Figure 2 
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Figure 3 
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Figure 4 
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Figure 5 
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Figure 6 
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Figure 7 
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Figure 8 
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Figure 9 
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Figure 10 
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EST OF RED VEL MADE AT TIME 6 
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Figure 11 
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Figure 12 
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Figure 13 
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Figure 14 
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Figure 15 
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Figure 16 
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Figure 18 
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Figure 22 
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Figure 23 
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Figure 25 
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Figure 26 
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Figure 27 
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APPENDIX A 

Node Occupancy and Its Perception, and COA Inference 

Describe the occupancy of a node, denoted by n, by the numbers of units of different 
types (e.g., heavy armor brigades, light armor brigades, etc.); £//(«; t ) is the number of 
such units of type i (i = 0, 1,2, I) at time t at node n; Uo(n; t ) can designate the empty 
node if necessary. 

Within a unit of type i, there may be several asset types. Agree that distinguishable 

assets occur in J classes, and that a unit of type i has a mean number of assets of type j 
(/= 1, 2, ..., J) equal to ctij(n,t), and a variance ajj(a;n,t) . Furthermore, adopt the 

provisional model that the actual number of assets of type j possessed by a particular 
randomly selected unit is a random variable, A(i,j, n, t ) with distribution function 
Fjj(x; n, t ). When convenient, and for illustration, we take n, t ), k= 1,2,..., 17/, 

i.e., the numbers of assets of type j owned by the 17/ copies on Node n to be independent 
and normally/Gaussian distributed. The time-dependent parameters ay and can reflect 

the fact that a campaign has been in progress for some time and attrition has occurred and 
is subject to change. Let 

_ u ‘ 

A(i,j,n,t) = (A.l) 

denote the total number of assets of type j possessed by units of type i at the n th node at 
time t. 

It is convenient to refer to the vector of distributions of typical asset-type numbers for 
a particular unit type as the signature (or asset signature) of the unit type. Note that 
signatures of different individual units of the same type will inevitably differ if their 
various asset counts differ, as could well happen. Let 
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(A.2) 


_ i u t 

1=1 1 =1 

the total number of assets of type j at the node at time t. 

Finally, let Sj(s, n; t) denote the total count of assets of type j by sensors of type s 
(s = 1,2, ...,S) at node n at time t. Sj represents the quantitative perception of the 
opponent’s type j asset level. It is convenient, although possibly optimistic, to assume 
initially that 

n;f)|^(y,w,r)] = A(j,n,t) (A3) 

with a variance that can depend upon sensor type s, node identity, n, and the time 
(especially if at night). Other influences on the ability of a sensor to assess the presence 
and strength of an opponent’s assets are the latter’s activity: radiation using radars, or 
communication attempts are two such clue-providers. 

Towards Probabilistic Node-Asset Assessment, Given Perception 

In an idealistic and theoretical sense it would be desirable for a protagonist in a 
theater-level campaign to possess the joint probability distribution of the types of units 
occupying all relevant nodes, along with their asset positions, at every time t. This 
distribution could be influenced by the protagonist’s disposition of its sensor assets and 
own force maneuver as well as by the opponent’s maneuver, presumably influenced by its 
own sensor asset disposition. Of course the raw sensor data can, and almost certainly will, 
be affected by deceptive tactics employed by the other side. Consequently fusion of the 
data and subsequent actions must attempt to wisely account for such a likelihood when 
situation assessments are made and future actions planned. 

The development of an automated systemic model of a theater campaign, or even a 
modest fragment thereof, must realistically stop short of even a completely integrated 
joint probabilistic assessment of nodal state for a particular protagonist. Instead, we will 
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endeavor to develop plausible marginal probabilistic descriptions of individual nodal 
states. 
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A Univariate Model 


We start with a univariate model, both for assets and sensors. For simplicity we drop 
explicit reference to n and t. Suppose that the node is occupied by (possibly multiple) 
units of one type, z, and a sensor of type 5 estimates the number of assets of type j at that 
node. We assume the following model: First, sensor assessment is randomly distributed 
around the actual asset number: 

.P{,Sy(.s;t) z.dx\A(j) = a} = —j== -expj-^x-zz) 2 / r 2 Syj \dx (A.4) 

V 2*7tT 


Second, asset number is randomly distributed around a known value, ayin, t), where the 
latter is derivable from the Table of Organization and Equipment (TOE), plus any 
information available concerning unit campaign experience (the unit may have been 
attrited). Assuming that actual y'^-asset values are independently distributed around the 
TOE value, then 


P{A{j) eda\Ui = u,U k = 0,k* z} 




a-uccij 


) 2 / uajj{a^da (A.5) 


It can easily be shown that 

P{<Sy [s, t) € dxpi = u,U k = 0, k * *} = 


rexpi 


V2^/«o-? (a) + 4 2ua 2 (a)+r SJ 


* (x-uccj/) 2 \dx 


(A.6) 


Further, by Bayes’ formula, 

P[A{j) eda\Sj(s;t) = x,U t = u,TJ k = 0,k * z'j 

1 f 1 

= — r — - -- exp<! — 

-v/2 7tVj[s,i,u,x) [ 2 

where 


( a-mj(s,i,u,x)^ 2 / vj(s,i,u,x)ida 
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Vj (s,i,u,x) = 


U<T\ 


|(a)4 


uafj(a) + rl 


(A. 8) 


/ . \ uafj(a) 

Mj{S,l,U,X) = — YTTl X ' 


■uau . 


(A.9) 


2 / \ 2 2/ \ 2 !/ 
wo-y(flj + r^ wcTy^J + r^ 

It next becomes necessary to account for the unknown number, u, of /-type units assumed 
to be present in order to assess the asset count distribution. To do so, let 


n(z,w) = P{Uj = u,U k = 0 ,k * /} (A. 10) 

the prior probability that there are u units of just one type, i, at the node. Recall that we 
are assuming there is only one type of unit at the node; thus, 

ES n ('>H (aid 

1=1 u 


We must compute 

n(/,u;x) = P^Ui = u,U k = 0 ,k * i\Sj(s;t) = xj 

= | P[U t =u,U k = 0 ,k*i,Aj e ( 5 ; t ) = xj (A. 12) 

= cn(z>) 1 exp{-J 27 ~y ~~ ~2 {x-uaijf) 

V2^V M<7 ?( a ) + 4- ua^aj + Tsj 

where the constant c is determined by the normalization condition. 

XZ n 0>;*) = 1 (A-13) 


Finally, then, the distribution of assets of type j is a probability (or convex) 
combination of normals: 

p\A(])^da\S,(s,l) = x' r 

s 

i,u 

and the joint distribution 
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(A. 15) 


P{Ui = u,Uk = 0,k*i,A(j) eda|iS)(.s;*) = x} 

= n(/, u;x ) - r— 1 —rexp j -\(a-mj (5, i,u,x)f /v) (5, i,u, x) W 

V2 xd j(s,i,u,x) [ 2 1 J 

This latter joint distribution is of the same form as the joint distribution prior to the sensor 
estimate. Thus similar calculations can be used to update the distribution as new sensor 
observations come in. To be specific let 

n(z',w;x(r),7) = P^Ui =u,U k = 0 ,k * = x(w),r(w) < (A.16) 

be the conditional distribution of the type of unit and number of that type of unit present 

given all the sensor observations that have occurred before time t. Let 
mj(t) = nij(s,i,u,x;t) (respectively o)(t) = v 2 j{s,i,u,x,t)) be the conditional expected 

value (respectively variance) of the number of asset j given there are u units of type i 
present at the node and all of the sensor observations before time t. Suppose a new sensor 
observation x(t + 1) from a sensor of type 5 occurs at t + 1. The posterior distribution n 
is updated as follows 

n(z',M;jc(r + l),r + l) 

/ 05 x (A. 17) 

= Cri(i,u;x(t),t)^x(t + l);m J (s > i,u,x(t);t),^u 2 j (s,i > u,x(t);t)+ rl) ' J 

where ^x; m, v ) is the normal density function with mean m and standard deviation v 
evaluated at x and C is a normalizing constant. Further the conditional moments of the 
number of asset j given the sensor observations and the type and number of unit present 
are updated as follows 
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nij(t + \) 


X (‘ + 1)+ -T77r. 72 m A>) 

U J 

and 


v){t)+T 2 sj " v *' v){t) + r]j 


(A. 18) 


Oj[t) + T s j 

An Independent Multivariate Model 

Once again assume that the node is occupied by (possibly multiple) units of one type 
and a sensor of type 5 estimates the number of assets of various types (say 2) at the node. 
Consider the following model. First a sensor assesses the numbers of different assets 
independently and each sensor asset assessment is randomly distributed around the actual 
asset number. 


(A. 19) 


P{5i(s;f) €<fri, 52 ( 3 *) edx 2 \A(l) = auA(2) = a 2 } 

2 j f 2 1 

=U-pr — exp r2 ( x j~ a j ) 7 

j =1 V27TT 5 j 

Second, given the type of unit and the number of units present, the counts of the different 
types of assets are independently distributed around the TOE values. 

P{^4(l) eda u A(2) eda 2 \Ui =u,Uk =0,fc*i} 

2 1 2 , (A.20) 

= n r 1 ~r — 7 ~\ ex Pl~ 2 ( a J - ua ») 

-42n<uc>ij(a) 1 J 

One can calculate. 


P{Sy(s;f) edxj-j = 1,2 |Ui = u,U k = 0,k * z'} 

= \\P{Sj{s\t)edXj-,A(j)edaj\j = \,2\Ui=u,U k =0,k*i} (A.21) 

a\a 2 
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Thus, if 


II (i,u) = P{U t =u,U k =0,k * i} 


(A.22) 


as before, then 


P^Ui =u,Uk =0,k * zjSi ( s ; t) = x\,S 2 ( s;t) = xi j 



= Tl{i,u\x\,X2). 


(A.23) 


Further 


P{A(j) edaj-J = 1,2|C7, = u,U k = 0 ,k * i,Sj(s,t) = xyj = 1,2} 


x 1 2 21 

= n nr~ 7 --exp{-}(a 7 -m y (5,z,w,xy)) /vj(s,i,u,xj) \daj 

j= i ■d2xVj(s,i,u,Xj) J 


(A.24) 


where 


Letting 


vj(s,i,u,xj) = 


ua}(a)T# 

uaf(a)+Tsj 


mj(s,i,u,Xj) = 


ua H a ) 4 

u<yjj{a) + t]j J u<jg(a) + r 2 sj 


UCCij 


(A.25) 



(A.26) 


the posterior distribution 
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p\a{J) edajij = 1,2, {7, = u,U k = 0 ,k * z|5\(.s;/) = xyj = 1,2} 

= Cn(i‘,w)]^| j;nij (s, i,u,xj),Vj(s,i,u, Xj ))^Xj ; uay , (ua] + r J )° 5 jdfly (A.27) 


= n(/, w; xi , x 2 £(a y ; mj ( s,i,u,xj ), vj(s,i,u,x y ))da y 

j= i 


where 


v 2 j(s,i,u,Xj) = 
mj(s,i,u,Xj) = 


wcri 




wcr-(a) + 4' 
ucrjj(a) 


Xt +■ 


ucrfj (a) + t]j u a 2 (a) + r 


2 ■««{/• 


(A.28) 


Recall that, 


= «,*7* =0 ,k*i,Aj edaj‘J- 1,2} 

2 (A.29) 

= n(i, y; M (Xjj , u ajj (a tya y 

y=i 

and thus the posterior distribution is of the same form as the prior. Hence, one can use the 
above posterior distribution as the new prior and compute a new posterior distribution 
when new sensor information becomes available. To be specific, let 
Ti(i,u\x{t),t) 

. (A. 30) 

= P{Uj = u,U k = 0 ,k* z|Sy(s(w),r(w)) = x(w),/(w) < t,j = 1,2} 

be the conditional distribution of the type of unit and the number of that type of unit 

present at the node given all the sensor observations that have occurred before time t. Let 
mj(t) = nij(s,i,u,x;t) (respectively v 2 (t) = Vj(s,i,u,x,t) ) be the conditional expected 

value (respectively variance) of the number of asset j given there are u units of type i 
present at the node and all of the sensor observations before time t. Suppose a new sensor 
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observation {xj{t +1)} from a sensor of type s occurs at t + 1. The posterior distribution 
n is updated as follows 

n(z,M;jc(t + l),t + l) 

where the product is over those types of assets for which there is an 
observation; m, v) is the normal density function with mean m and standard deviation 
v evaluated at x and C is a normalizing constant. Further the conditional moments of the 
number of asset j given the sensor observations and the type and number of unit present 
are updated as follows 


= Cn(i,u\x{t),t)Y\^Xj{t + ^),mj{sa,u,x\t), v 2 {t) + T 2 S j ° 


(A.31) 



v 2 (t + l) = 


v 2 j{t) + r 2 sj 


Example 

In this example there are 2 types of units: I and II. There are 3 asset types. Suppose 
unit I has 0 assets of type 1, 100 assets of type 2 and 30 assets of type 3. Suppose unit II 
has 50 assets of type 1,150 assets of type 2 and 50 assets of type 3. 

Suppose there are 2 sensor types. Sensor type 1 can estimate the count of assets of 
type 1 and 2, but not 3. Sensor type 3 can estimate the count of assets of type 3 but not 1 
and 2. 

Assume the following prior distribution. 


113 


where 


P{ A {j) eda 

,\Ui = 

u;U k 

= 0, k ^ zj = 

= 



n-jL 1 

7=7 V 2n 'Ju<j i j(a) 

exp{- 

2 ( a j ~ ua v 

) 2 /W<7,}(a)jj<2y 

P{Ui 

= u,U k 

= 0,k 

^‘} = 

= — for i 
4 

= 1,2; 

u = 1,2 



0 

if j 

= 1 


0.5 

if j 

= 1 

«i j = < 

100 

if j 

= 2 


25 

1 

if j 

= 2 


30 

if j 

= 3 


l 

5 

if j 

= 3 


'75 

if j 

= 1 


'10 

if j 

= 1 

a 2J =< 

150 

if j 

= 2 

<*2 j=< 

25 

if j 

= 2 


50 

if j 

= 3 


5 

if j 

= 3 


(A.33) 


(A-3 4) 


that is, the number of assets of each type are independently normally distributed and the 
prior distribution of the unit type and number of units occupying the node assigns 
probability 0.25 to each of the events: there is one (respectively two) units of type I; there 
is one (respectively two) units of type H 

Assume the distribution of the sensor estimate of the count of assets at the node is 
P{Sy(l,n;f) edxj-J = 1,2| A(j) = a/J = 1,2} 

7sR exp B 

where 



(xj-cij) 2 1 rlj^dxj 






'25 

35 


if y = l 
if 7 = 2 


P{5 3 (2, n,t) e dx 2 \ A 3 = a 3 } = exp j- ^(*3 - a 3 ) / t 2 3 1 


(A. 3 6) 
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where ^3 = 10 . 

The results of the updating procedure appears in Table A1. Displayed are the time of 
the sensor report, the type of sensor making the report, the estimated number of assets of 
each type made by the sensor, the posterior distribution of the type of unit and number of 
units present, and the conditional posterior means and standard deviations of the number 
of assets of each type present given the type of unit and the number of units present. 

For time 0, the prior distribution of the unit type and number of units present 
ri(/, u) = 1/4 is displayed. Also displayed are the conditional moments of the prior normal 
distribution of numbers of assets of each type given the unit type and number of units 
present. At time 1 there is a report from a sensor of type 1 that there are 10 assets of 
type 1 and 80 assets of type 2. This information updates the probability that there is 1 
unit of type I at the node to 0.94. At time 2 there is a report from a sensor of type 2 that 
there are 35 assets of type 3. This information updates the probability that there is 1 unit 
of type I at the node to 0.99. At time 3 there is a report from a sensor of type 1 that there 
are 60 assets of type 1 and 140 assets of type 2 . This information decreases the probability 
that there is one unit of type I at the node to 0.81 and increases the probability that there is 
one unit of type 2 at the node to 0.18. At time 4 there is a report from a sensor of type 1 
that there are 90 assets of type 1 and 200 assets of type 2. This information increases the 
probability that there is one unit of type II at the node to 1. At time 5 there is a report 
from a sensor of type 2 that there are 70 assets of type 3. This information decreases the 
probability that there is one unit of type II at the node to 0.13 and increases the probability 
that there are 2 units of type II at the node to 0.87. 

The sensor observations were chosen so that the observations at times 1 and 2 mi ght 
come from one unit of type I. The sensor observations after time 2 may come from 
another situation. Note that the procedure does respond to the change. 
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To further investigate the responsiveness of the procedure to change, a simulation was 
conducted. The ground truth for the simulation is as follows. For the first 5 time periods 
there is 1 unit of type I at the node. The (true) count of each asset type is drawn from the 
prior normal distribution for one unit of type I and rounded down to the closest integer; if 
the simulated asset count is negative it is set equal to zero. The (true) count of asset 1 
(respectively 2, 3) is 0 (respectively 69, 38) for the first five time periods. For the 
remainder of the time periods, there are 2 units of type II at the node. Once again the 
count of assets of each type is drawn from the prior normal distribution for two units of 
type II and rounded down to the closest integer. The resulting count of asset 1 
(respectively 2, 3) is 149 (respectively 290, 108). The standard deviation of sensor 1 on 
asset 1 (respectively 2) is 25 (respectively 35). The standard deviation of sensor 2 on asset 
3 is 10. 

For each time period, a random number is drawn to determine the type of sensor 
making the observation; (each sensor type is equally likely). Given the sensor type, the 
sensor observation is drawn from a normal distribution with mean the true asset number 
and standard deviation for the asset type and sensor type; the resulting number is set equal 
to zero if it is negative. 

The results of the simulation appear in Table A2. Once again the table lists the times 
of the observations, the type of sensor making the observation, the sensor’s estimated 
number of assets, the posterior distribution of the unit type and number of units after each 
observation, and the conditional posterior moments. The row for time 0 displays the prior 
distribution. 

Recall that there is one unit of type I at the node for times 1-5. The posterior 
distributions quickly reflect this ground truth; after the first observation the posterior 
probability is 0.89 that there is one unit of type I; after the second observation, this 
posterior probability is increased to 0.99; after the third observation to 1.0. 
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Recall that for times 6 on, there are two units of type II at the node. The observation at 
time 6 causes the posterior distribution of there being 1 unit of type I to decrease to 0.05 
and increases the probability of there being 2 units of type I to 0.95. Thus, the updating 
procedure does reflect the fact that a change has occurred but does not correctly identify 
the type and number of units. It is not until after the observation at time 12, that the 
updating procedure puts a sizable posterior probability on the event that there are 2 units 
of type II. This sluggish behavior is probability due to that fact that after the observation 
at time 5, the posterior probability that there are two units of type II at the node is (a very 
unlikely) 10* 26 . Subsequent observations raise the probability to 10' 5 after the 11 th 
observation. Thus, it appears that the Bayesian methodology must overcome its belief 
after time 5 that it is extremely unlikely that there are 2 units of type II at the node. 

Some Generalizations 

Let X(n;t) = (Ui(n;t), U 2 (n;t) ..., Uj(n;ty) be a vector describing the number of units of 
each type at node n at time t; Ui(n,t ) is the number of units of type i at node n at time t. If 
the node is unoccupied at time t, X(n,t ) = (0,0,.. .,0). 

The total number of assets of type j at the node at time t is 

_ / U,(n,t) 

A(j,n,t) = '£ {i,j,n,t). 

i=l k =1 

Assume {A]^ij,n,t)} are independent normal random variables with mean ay and 
variance crfj(a ). If Ui(n,t) > 0 for some i then A(j,n,t ) is normal with mean 

I Ui 

m j {ux...,ur,n$) = Y J Y J a ij 

i=i k=\ 

and variance 


1 Ui 

Vj(u x ...,u I \n\0) = Y J Yj °i ( a ) 


i=i *=i 
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If the node is empty we assume A[j,n,t ) is normal with mean 0 and small standard 
deviation, e.g., 0.25. 

Let 7iu\, M 2 , ...,m/, 0) denote the prior possibility that there are n z units of type i 
present at the node. 

Assume that a sensor of type 5 is able to estimate the number of assets of type j for 
jsJ(s) where J(s) is a subset of the assets types. For each sensor type s 


p{Sj(s-,n;t) edxj'J zJ(s\A(j,n,t) = a/,j e. J(j)|} 


-n 

jeJ(s) 








that is the sensor observations are independent normally distributed about the true number 
of asset j. 

Once again a recursive procedure for updating the distribution of the number and 
types of units at the node can be obtained; as well as a procedure for updating the 
conditional mean and variance of the number of asset j at the node given the number and 
types of units. The development is similar to that of the previous sections. To be specific 
let 


n((Mi,...,M/);x(r),r) 

= P{Ui= Mi,...,£7/ = M 7 |5 y (5(w),r(w)) = x(w),y ei(s),/(w)<t} 

be the conditional distribution of the type of unit and the number of that type of unit 

present at the node given all the sensor observations that have occurred before time t. Let 
mj((u i, t) (respectively o 2 j{u\,...,ui)\t ) be the conditional expected value 

(respectively variance) of the number of assets of typey given U\ =u \,..., Uj= uj and all 
the sensor observations before time t. Suppose a new sensor observation 
{xj(t + 1 );jeJ(s)} becomes available. The updated probability of unit type and number is 
computed as 
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n((«i ,..., M/); x(t), x(t + 1); t +1) 


0.5 


= Cn((«i,...,ii/);*(/);f) Yl QXj{t + l)\mj(t),(y 2 j(t)+T 2 sj ) 

jeJ(s) ' 

where C is a normalizing constant. The conditional moments are updated as follows. 
mj((ui,...,ui);t + 1) 


T v 


uj({ Ul 0 , _ 

v)[(u x ,...,u I )\t) + z 2 sj J wy((i/i,...,M/);r)+T; 




and 


^y'((wiv>M/);? + l) =-—- L - Z - 




+ r; 


A Modified Bayesian Approach 

The simulation results of Table A2 suggest that the complete Bayesian approach can 
be rather sluggish in identifying the true number and types of units at a node. Simulation 
experiments not reported here suggest that in some cases the procedure may never 
identify the true number and types of units. As mentioned earlier, one cause for this 
behavior appears to be that the posterior distribution of the number and types of units at 
the node can attribute essentially 0 probability to particular configurations; this 0 
probability can make it difficult for the procedure to adjust to new circumstances. 

Another aspect of the Bayesian procedure is that it continually updates the conditional 
means and variances of the amount of each type of asset for each unit type and number 
configuration. This continual updating may result in the posterior conditional moments 
being far away from the original conditional moments; this may negatively influence the 
ability of the Bayesian procedure to correctly identify the numbers and types of units at a 
node in a new situation. 

Preliminary experiments have suggested that the following may result in an improved 
ability to correctly identify the number and type of units at a node. At each time t compare 
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the just computed posterior distribution n((u\, ..., uj)\t) to the prior distribution 
<(tq, uj)\t-\). If 

max Jn(( Wl )jt)—n((wj.. ,M/); r —1)|>0.8, 

say, declare that a change has occurred. Use an equally likely prior distribution for unit 
types and numbers when the next sensor observation arrives; also reset the conditional 
moments to their original values before updating. 

An Example 

Tables A3 and A4 report the results of a simulation study. The scenario is as follows. 
One unit of type I has a mean number of asset 1 (respectively assets 2 and 3) of 0 
(respectively 100, 30) with standard deviation for the number of asset 1 (respectively 
assets 2 and 3) of 0.5 (respectively 25 and 5). One unit of type II has a mean number of 
asset 1 (respectively assets 2 and 3) of 75 (respectively 15 and 50) with standard 
deviation for the number of asset 1 (respectively assets 2 and 3) of 10 (respectively 25 
and 5). A sensor of type 1 can estimate a number of assets of type 1 and type 2 but not 
type 3; the error standard deviation is 25, (respectively 35) for type 1 (respectively type 2) 
assets. A sensor of type 2 can estimate numbers of assets of type 3 only with error 
standard deviation 10. 

The possibilities for node occupancy are (0, 0) (empty); (1, 0) (1 unit of type I); (2, 0) 
(2 units of type I); (0,1) (1 unit of type II); (0, 2) (2 units of type II), and (1,1) (1 unit of 
type I and 1 unit of type II). The prior distribution is that each of these configurations is 
equally likely. 

During the first 5 time periods there is 1 unit of type I at the node; the true number of 
assets of type 1 (respectively asset type 2 and 3) is 0 (respectively 69 and 39); these 
numbers are obtained by drawing from the appropriate normal distribution for asset type, 
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making the resulting number 0 if it is negative and rounding down to the next integer 
value. 

During times 6-15 there are 1 unit of type I and 1 unit of type II at the node with 99 
combined assets of type 1, 231 combined assets of type 2, and 63 combined assets of type 
3; these true numbers of assets are obtained in a similar manner to those of the first 5 time 
periods. 

At each time a sensor type is randomly selected and the sensor observation(s) are 
drawn from a normal distribution with mean the true number of assets of that type and 
standard deviation for that asset type for that sensor; the observations are set equal to 0 if 
they are negative. 

Table A3 displays results of using the Bayesian updating procedure. Displayed are the 
data, the posterior probabilities of each configuration, and the conditional mean number 
of each asset type given the configuration and the observations. The row at time 0 
displays the prior information. The posterior distributions indicate that the procedure uses 
the sensor observations to correctly classify by time 5 that there is 1 unit of type I at the 
node. By time 7 the posterior distribution indicates that there is evidence that the 
occupancy of the node has changed. However, the posterior distribution is never able to 
assign an appreciable probability to the correct identity of 1 unit of type I and 1 unit of 
type II. 

Table A4 displays the results of using the Bayesian updating procedure with the 
following adjustment. At each time the prior and posterior distribution are compared. If 
the total variation is less than 0.8, nothing is done. If it is greater than 0.8, the posterior is 
replaced with the original prior and the posterior conditional moments are replaced with 
the original prior conditional moments for the calculations the next time period. In Table 
A4, the total variation of the prior and posterior exceeded 0.8 at time 11. Thus the results 
for times 12-15 differ from those of Table A3. Note that in this case the posterior 
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distribution is able to assign a largish probability to the true configuration of 1 unit of 
type I and 1 unit of type n. 

Note that the random assignment of sensor type has assigned a type 2 sensor to make 
most of the observations. Recall that a type 2 sensor can only estimate numbers of assets 
of type 3. However, the most obvious feature of a type I unit is that it has no assets of 
type I. Thus, it is not surprising that the Bayesian procedure is not perfect. 

Note also, that even though the occupancy of the node changed at time 6, it was not 
until time 11 until the total variation between the prior and posterior indicates it; note that 
time 11 is also the first observation from a sensor of type I after time 5. 
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Table A1 
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Table A3 

Updating Perceptions at a Node; No Adjustment 
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Table A4 

Updating Perceptions at a Node; With Adjustment 
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